Changeset 357 for LMDZ.3.3/trunk/libf


Ignore:
Timestamp:
May 7, 2002, 3:04:11 PM (23 years ago)
Author:
lmdz
Message:

Remplacement, dans certaines circonstances l'ancienne version pouvait faire
planter le modele. FH
LF

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r2 r357  
    2626c   ----------
    2727      REAL masse(ip1jmp1,llm),pente_max
     28c      REAL masse(iip1,jjp1,llm),pente_max
    2829      REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
    2930      REAL q(ip1jmp1,llm)
     31c      REAL q(iip1,jjp1,llm)
    3032      REAL w(ip1jmp1,llm),pdt
    3133c
    … …  
    3335c   ---------
    3436c
    35       INTEGER ij,l
    36 c
    37       REAL zm(ip1jmp1,llm)
     37      INTEGER i,ij,l,j,ii
     38      INTEGER ijlqmin,iqmin,jqmin,lqmin
     39c
     40      REAL zm(ip1jmp1,llm),newmasse
    3841      REAL mu(ip1jmp1,llm)
    3942      REAL mv(ip1jm,llm)
    4043      REAL mw(ip1jmp1,llm+1)
    41       REAL zq(ip1jmp1,llm)
    42       REAL temps1,temps2,temps3
     44      REAL zq(ip1jmp1,llm),zz
     45      REAL dqx(ip1jmp1,llm),dqy(ip1jmp1,llm),dqz(ip1jmp1,llm)
     46      REAL second,temps0,temps1,temps2,temps3
     47      REAL ztemps1,ztemps2,ztemps3
    4348      REAL zzpbar, zzw
    4449      LOGICAL testcpu
    4550      SAVE testcpu
    4651      SAVE temps1,temps2,temps3
     52      INTEGER iminn,imaxx
    4753
    4854      REAL qmin,qmax
    … …  
    5157      DATA temps1,temps2,temps3/0.,0.,0./
    5258
    53 c      PRINT*,'Debut vlsplt version debug sans vly'
    5459
    5560        zzpbar = 0.5 * pdt
    … …  
    7075         mw(ij,llm+1)=0.
    7176      ENDDO
    72 
     77     
    7378      CALL SCOPY(ijp1llm,q,1,zq,1)
    7479      CALL SCOPY(ijp1llm,masse,1,zm,1)
    7580
    76 c      call minmaxq(zq,qmin,qmax,'avant vlx     ')
     81        print*,'Entree vlx1'
     82c       call minmaxq(zq,qmin,qmax,'avant vlx     ')
    7783      call vlx(zq,pente_max,zm,mu)
    78 
    79 
    80 c     call minmaxq(zq,qmin,qmax,'avant vly     ')
    81 
     84        print*,'Sortie vlx1'
     85c       call minmaxq(zq,qmin,qmax,'apres vlx1    ')
     86
     87         print*,'Entree vly1'
    8288      call vly(zq,pente_max,zm,mv)
    83 
    84 
    85 c      call minmaxq(zq,qmin,qmax,'avant vlz     ')
    86 
     89c       call minmaxq(zq,qmin,qmax,'apres vly1     ')
     90        print*,'Sortie vly1'
    8791      call vlz(zq,pente_max,zm,mw)
    88 
    89 
    90 c     call minmaxq(zq,qmin,qmax,'avant vly     ')
    91 c     call minmaxq(zm,qmin,qmax,'M avant vly     ')
     92c       call minmaxq(zq,qmin,qmax,'apres vlz     ')
     93
    9294
    9395      call vly(zq,pente_max,zm,mv)
    94 
    95 
    96 c     call minmaxq(zq,qmin,qmax,'avant vlx     ')
    97 c     call minmaxq(zm,qmin,qmax,'M avant vlx     ')
     96c       call minmaxq(zq,qmin,qmax,'apres vly     ')
     97
    9898
    9999      call vlx(zq,pente_max,zm,mu)
    100 
    101 c     call minmaxq(zq,qmin,qmax,'apres vlx     ')
    102 c     call minmaxq(zm,qmin,qmax,'M apres vlx     ')
    103 
     100c       call minmaxq(zq,qmin,qmax,'apres vlx2    ')
     101       
    104102
    105103      DO l=1,llm
    … …  
    115113      END
    116114      SUBROUTINE vlx(q,pente_max,masse,u_m)
    117 c
     115
    118116c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
    119117c
    … …  
    137135c   ----------
    138136      REAL masse(ip1jmp1,llm),pente_max
    139       REAL u_m( ip1jmp1,llm )
     137      REAL u_m( ip1jmp1,llm ),pbarv( iip1,jjm,llm)
    140138      REAL q(ip1jmp1,llm)
     139      REAL w(ip1jmp1,llm)
    141140c
    142141c      Local
    … …  
    147146c
    148147      REAL new_m,zu_m,zdum(ip1jmp1,llm)
    149       REAL dxq(ip1jmp1,llm),dxqu(ip1jmp1)
     148      REAL sigu(ip1jmp1),dxq(ip1jmp1,llm),dxqu(ip1jmp1)
    150149      REAL zz(ip1jmp1)
    151150      REAL adxqu(ip1jmp1),dxqmax(ip1jmp1,llm)
    152151      REAL u_mq(ip1jmp1,llm)
    153152
    154       Logical first,testcpu
     153      Logical extremum,first,testcpu
    155154      SAVE first,testcpu
    156155
    157156      REAL      SSUM
    158157      EXTERNAL  SSUM
    159       REAL temps1,temps2,temps3,temps4,temps5
    160       SAVE temps1,temps2,temps3,temps4,temps5
    161 
     158      REAL temps0,temps1,temps2,temps3,temps4,temps5,second
     159      SAVE temps0,temps1,temps2,temps3,temps4,temps5
     160
     161      REAL z1,z2,z3
    162162
    163163      DATA first,testcpu/.true.,.false./
    … …  
    176176
    177177      IF (pente_max.gt.-1.e-5) THEN
    178 c     IF (pente_max.gt.10) THEN
     178c       IF (pente_max.gt.10) THEN
    179179
    180180c   calcul des pentes avec limitation, Van Leer scheme I:
    … …  
    230230
    231231         ENDDO ! l=1,llm
     232        print*,'Ok calcul des pentes'
    232233
    233234      ELSE ! (pente_max.lt.-1.e-5)
    … …  
    266267            dxq(ij-iim,l)=dxq(ij,l)
    267268         ENDDO
    268 
    269269         DO ij=1,ip1jmp1
    270270            iadvplus(ij,l)=0
    … …  
    273273      ENDDO
    274274
     275         print*,'Bouclage en iip1'
    275276
    276277c   calcul des flux a gauche et a droite
    … …  
    294295c   on cumule le flux correspondant a toutes les mailles dont la masse
    295296c   au travers de la paroi pENDant le pas de temps.
     297        print*,'Cumule ....'
     298
    296299      DO l=1,llm
    297300       DO ij=iip2,ip1jm-1
     301c       print*,'masse(',ij,')=',masse(ij,l)
    298302          IF (u_m(ij,l).gt.0.) THEN
    299303             zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l)
    … …  
    306310      ENDDO
    307311#endif
    308 
    309 
     312c       stop
     313
     314c       go to 9999
    310315c   detection des points ou on advecte plus que la masse de la
    311316c   maille
    … …  
    318323         ENDDO
    319324      ENDDO
     325        print*,'Ok test 1'
    320326      DO l=1,llm
    321327       DO ij=iip1+iip1,ip1jm,iip1
    … …  
    323329       ENDDO
    324330      ENDDO
    325 
     331         print*,'Ok test 2'
    326332
    327333
    … …  
    342348
    343349      IF(n0.gt.1) THEN
    344 ccc      PRINT*,'Nombre de points pour lesquels on advect plus que le'
    345 ccc     &       ,'contenu de la maille : ',n0
     350      PRINT*,'Nombre de points pour lesquels on advect plus que le'
     351     &       ,'contenu de la maille : ',n0
    346352
    347353         DO l=1,llm
    … …  
    395401         ENDDO
    396402      ENDIF  ! n0.gt.0
    397 
     4039999    continue
    398404
    399405
    400406c   bouclage en latitude
    401 
     407        print*,'Avant bouclage en latitude'
    402408      DO l=1,llm
    403409        DO ij=iip1+iip1,ip1jm,iip1
    … …  
    423429         ENDDO
    424430      ENDDO
    425 
    426431c     CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
    427432c     CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1)
    … …  
    456461      REAL masse(ip1jmp1,llm),pente_max
    457462      REAL masse_adv_v( ip1jm,llm)
    458       REAL q(ip1jmp1,llm)
     463      REAL q(ip1jmp1,llm), dq( ip1jmp1,llm)
    459464c
    460465c      Local
    … …  
    464469c
    465470      REAL airej2,airejjm,airescb(iim),airesch(iim)
    466       REAL dyq(ip1jmp1,llm),dyqv(ip1jm)
     471      REAL dyq(ip1jmp1,llm),dyqv(ip1jm),zdvm(ip1jmp1,llm)
    467472      REAL adyqv(ip1jm),dyqmax(ip1jmp1)
    468473      REAL qbyv(ip1jm,llm)
    469474
    470       REAL qpns,qpsn,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
     475      REAL qpns,qpsn,apn,aps,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
    471476c     REAL newq,oldmasse
    472       Logical first,testcpu
    473       REAL temps0,temps1,temps2,temps3,temps4,temps5
     477      Logical extremum,first,testcpu
     478      REAL temps0,temps1,temps2,temps3,temps4,temps5,second
    474479      SAVE temps0,temps1,temps2,temps3,temps4,temps5
    475480      SAVE first,testcpu
    476481
    477482      REAL convpn,convps,convmpn,convmps
     483      real massepn,masseps,qpn,qps
    478484      REAL sinlon(iip1),sinlondlon(iip1)
    479485      REAL coslon(iip1),coslondlon(iip1)
    … …  
    506512
    507513c
    508 
     514        PRINT*,'CALCUL EN LATITUDE'
    509515
    510516      DO l = 1, llm
    … …  
    565571c   calcul des pentes limites aux poles
    566572
     573      goto 8888
    567574      fn=1.
    568575      fs=1.
    … …  
    578585         dyq(ij,l)=fn*dyq(ij,l)
    579586         dyq(ip1jm+ij,l)=fs*dyq(ip1jm+ij,l)
     587      ENDDO
     5888888    continue
     589      DO ij=1,iip1
     590         dyq(ij,l)=0.
     591         dyq(ip1jm+ij,l)=0.
    580592      ENDDO
    581593
    … …  
    682694         ENDDO
    683695c.-. ancienne version
    684          convpn=SSUM(iim,qbyv(1,l),1)/apoln
    685          convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln
    686          DO ij = 1,iip1
    687             newmasse=masse(ij,l)+convmpn*aire(ij)
    688             q(ij,l)=(q(ij,l)*masse(ij,l)+convpn*aire(ij))/
    689      &               newmasse
    690             masse(ij,l)=newmasse
    691          ENDDO
    692          convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)/apols
    693          convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)/apols
    694          DO ij = ip1jm+1,ip1jmp1
    695             newmasse=masse(ij,l)+convmps*aire(ij)
    696             q(ij,l)=(q(ij,l)*masse(ij,l)+convps*aire(ij))/
    697      &               newmasse
    698             masse(ij,l)=newmasse
    699          ENDDO
     696c        convpn=SSUM(iim,qbyv(1,l),1)/apoln
     697c        convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln
     698
     699         convpn=SSUM(iim,qbyv(1,l),1)
     700         convmpn=ssum(iim,masse_adv_v(1,l),1)
     701         massepn=ssum(iim,masse(1,l),1)
     702         qpn=0.
     703         do ij=1,iim
     704            qpn=qpn+masse(ij,l)*q(ij,l)
     705         enddo
     706         qpn=(qpn+convpn)/(massepn+convmpn)
     707         do ij=1,iip1
     708            q(ij,l)=qpn
     709         enddo
     710
     711c        convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)/apols
     712c        convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)/apols
     713
     714         convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
     715         convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
     716         masseps=ssum(iim, masse(ip1jm+1,l),1)
     717         qps=0.
     718         do ij = ip1jm+1,ip1jmp1-1
     719            qps=qps+masse(ij,l)*q(ij,l)
     720         enddo
     721         qps=(qps+convps)/(masseps+convmps)
     722         do ij=ip1jm+1,ip1jmp1
     723            q(ij,l)=qps
     724         enddo
     725
    700726c.-. fin ancienne version
    701727
    … …  
    852878      RETURN
    853879      END
    854       SUBROUTINE minmaxq(zq,qmin,qmax,comment)
     880c      SUBROUTINE minmaxq(zq,qmin,qmax,comment)
     881c
     882c#include "dimensions.h"
     883c#include "paramet.h"
     884
     885c      CHARACTER*(*) comment
     886c      real qmin,qmax
     887c      real zq(ip1jmp1,llm)
     888
     889c      INTEGER jadrs(ip1jmp1), jbad, k, i
     890
     891
     892c      DO k = 1, llm
     893c         jbad = 0
     894c         DO i = 1, ip1jmp1
     895c         IF (zq(i,k).GT.qmax .OR. zq(i,k).LT.qmin) THEN
     896c            jbad = jbad + 1
     897c            jadrs(jbad) = i
     898c         ENDIF
     899c         ENDDO
     900c         IF (jbad.GT.0) THEN
     901c         PRINT*, comment
     902c         DO i = 1, jbad
     903cc            PRINT*, "i,k,zq=", jadrs(i),k,zq(jadrs(i),k)
     904c         ENDDO
     905c         ENDIF
     906c      ENDDO
     907
     908c      return
     909c      end
     910      subroutine minmaxq(zq,qmin,qmax,comment)
    855911
    856912#include "dimensions.h"
    857913#include "paramet.h"
    858914
    859       CHARACTER*(*) comment
     915      character*20 comment
    860916      real qmin,qmax
    861917      real zq(ip1jmp1,llm)
    862 
    863       INTEGER jadrs(ip1jmp1), jbad, k, i
    864 
    865       DO k = 1, llm
    866          jbad = 0
    867          DO i = 1, ip1jmp1
    868          IF (zq(i,k).GT.qmax .OR. zq(i,k).LT.qmin) THEN
    869             jbad = jbad + 1
    870             jadrs(jbad) = i
    871          ENDIF
    872          ENDDO
    873          IF (jbad.GT.0) THEN
    874          PRINT*, comment
    875          DO i = 1, jbad
    876 cc            PRINT*, "i,k,zq=", jadrs(i),k,zq(jadrs(i),k)
    877          ENDDO
    878          ENDIF
    879       ENDDO
     918      real zzq(iip1,jjp1,llm)
     919
     920      integer imin,jmin,lmin,ijlmin
     921      integer imax,jmax,lmax,ijlmax
     922
     923      integer ismin,ismax
     924
     925      call scopy (ip1jmp1*llm,zq,1,zzq,1)
     926
     927      ijlmin=ismin(ijp1llm,zq,1)
     928      lmin=(ijlmin-1)/ip1jmp1+1
     929      ijlmin=ijlmin-(lmin-1.)*ip1jmp1
     930      jmin=(ijlmin-1)/iip1+1
     931      imin=ijlmin-(jmin-1.)*iip1
     932      zqmin=zq(ijlmin,lmin)
     933
     934      ijlmax=ismax(ijp1llm,zq,1)
     935      lmax=(ijlmax-1)/ip1jmp1+1
     936      ijlmax=ijlmax-(lmax-1.)*ip1jmp1
     937      jmax=(ijlmax-1)/iip1+1
     938      imax=ijlmax-(jmax-1.)*iip1
     939      zqmax=zq(ijlmax,lmax)
     940
     941       if(zqmin.lt.qmin)
     942c     s     write(*,9999) comment,
     943     s     write(*,*) comment,
     944     s     imin,jmin,lmin,zqmin,zzq(imin,jmin,lmin)
     945       if(zqmax.gt.qmax)
     946c     s     write(*,9999) comment,
     947     s     write(*,*) comment,
     948     s     imax,jmax,lmax,zqmax,zzq(imax,jmax,lmax)
    880949
    881950      return
     9519999  format(a20,'  q(',i3,',',i2,',',i2,')=',e12.5,e12.5)
    882952      end
    883953
Note: See TracChangeset for help on using the changeset viewer.