/[lmdze]/trunk/Sources/dyn3d/PPM3d/ppm3d.f
ViewVC logotype

Diff of /trunk/Sources/dyn3d/PPM3d/ppm3d.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 10 by guez, Fri Apr 18 14:45:53 2008 UTC revision 32 by guez, Tue Apr 6 17:52:58 2010 UTC
# Line 397  C Line 397  C
397        do J=2,JMR        do J=2,JMR
398        DTDX(j)  = DT / ( DL*AE*COSP(J) )        DTDX(j)  = DT / ( DL*AE*COSP(J) )
399    
 c     print*,'dtdx=',dtdx(j)  
400        DTDX5(j) = 0.5*DTDX(j)        DTDX5(j) = 0.5*DTDX(j)
401        enddo        enddo
402  C  C
403                
404        DTDY  = DT /(AE*DP)        DTDY  = DT /(AE*DP)
 c      print*,'dtdy=',dtdy  
405        DTDY5 = 0.5*DTDY        DTDY5 = 0.5*DTDY
406  C  C
 c      write(6,*) 'J1=',J1,' J2=', J2  
407        endif        endif
408  C  C
409  C *********** End Initialization **********************  C *********** End Initialization **********************
# Line 726  C Line 723  C
723        DO j=1,JNP        DO j=1,JNP
724        DO i=1,IMR        DO i=1,IMR
725              Q(i,j,k,IC) = DQ(i,j,k,IC) / delp2(i,j,k)              Q(i,j,k,IC) = DQ(i,j,k,IC) / delp2(i,j,k)
 c            print*,'i=',i,'j=',j,'k=',k,'Q(i,j,k,IC)=',Q(i,j,k,IC)  
726        enddo        enddo
727        enddo        enddo
728        enddo        enddo
# Line 863  C Line 859  C
859        AL(i,k) = wk1(i,k-1) + c1 + c2 *        AL(i,k) = wk1(i,k-1) + c1 + c2 *
860       &        ( wk2(i,k  )*(c1*(A1 - A2)+A2*flux(i,k-1)) -       &        ( wk2(i,k  )*(c1*(A1 - A2)+A2*flux(i,k-1)) -
861       &          wk2(i,k-1)*A1*flux(i,k)  )       &          wk2(i,k-1)*A1*flux(i,k)  )
 C      print *,'AL1',i,k, AL(i,k)  
862  12    CONTINUE  12    CONTINUE
863  14    continue  14    continue
864  C  C
865        do 20 i=1,IMR*NLAYM1        do 20 i=1,IMR*NLAYM1
866        AR(i,1) = AL(i,2)        AR(i,1) = AL(i,2)
 C      print *,'AR1',i,AR(i,1)  
867  20    continue  20    continue
868  C  C
869        do 30 i=1,IMR*NLAY        do 30 i=1,IMR*NLAY
870        A6(i,1) = 3.*(wk1(i,1)+wk1(i,1) - (AL(i,1)+AR(i,1)))        A6(i,1) = 3.*(wk1(i,1)+wk1(i,1) - (AL(i,1)+AR(i,1)))
 C      print *,'A61',i,A6(i,1)  
871  30    continue  30    continue
872  C  C
873  C****6***0*********0*********0*********0*********0*********0**********72  C****6***0*********0*********0*********0*********0*********0**********72
# Line 895  C Line 888  C
888          CM = wz2(i,1) / wk2(i,1)          CM = wz2(i,1) / wk2(i,1)
889          flux(i,2) = AR(i,1)+0.5*CM*(AL(i,1)-AR(i,1)+A6(i,1)*(1.-R23*CM))          flux(i,2) = AR(i,1)+0.5*CM*(AL(i,1)-AR(i,1)+A6(i,1)*(1.-R23*CM))
890        else        else
 C        print *,'test2-0',i,j,wz2(i,1),wk2(i,2)  
891          CP= wz2(i,1) / wk2(i,2)                  CP= wz2(i,1) / wk2(i,2)        
 C        print *,'testCP',CP  
892          flux(i,2) = AL(i,2)+0.5*CP*(AL(i,2)-AR(i,2)-A6(i,2)*(1.+R23*CP))          flux(i,2) = AL(i,2)+0.5*CP*(AL(i,2)-AR(i,2)-A6(i,2)*(1.+R23*CP))
 C        print *,'test2',i, AL(i,2),AR(i,2),A6(i,2),R23  
893        endif        endif
894  140   continue  140   continue
895  C  C
# Line 1054  c            LMT = IORD - 3 Line 1044  c            LMT = IORD - 3
1044  c      endif  c      endif
1045  C  C
1046        LMT = IORD - 3        LMT = IORD - 3
1047  c      write(6,*) 'PPM option in E-W direction = ', LMT  
 c      first = .false.  
 C      endif  
 C  
1048        DO 10 i=1,IMR        DO 10 i=1,IMR
1049  10    AL(i) = 0.5*(p(i-1)+p(i)) + (DC(i-1) - DC(i))*R3  10    AL(i) = 0.5*(p(i-1)+p(i)) + (DC(i-1) - DC(i))*R3
1050  C  C
# Line 1335  C Poles: Line 1322  C Poles:
1322          wk(i,JNP+2) = p(i+IMH,JNP-2)          wk(i,JNP+2) = p(i+IMH,JNP-2)
1323          wk(i+IMH,JNP+2) = p(i,JNP-2)          wk(i+IMH,JNP+2) = p(i,JNP-2)
1324          enddo          enddo
1325  c        write(*,*) 'toto 1'  
 C --------------------------------  
1326        IF(IAD.eq.2) then        IF(IAD.eq.2) then
1327        do j=j1-1,j2+1        do j=j1-1,j2+1
1328        do i=1,IMR        do i=1,IMR
 c      write(*,*) 'avt NINT','i=',i,'j=',j  
1329        JP = NINT(VA(i,j))              JP = NINT(VA(i,j))      
1330        rv = JP - VA(i,j)        rv = JP - VA(i,j)
 c      write(*,*) 'VA=',VA(i,j), 'JP1=',JP,'rv=',rv  
1331        JP = j - JP        JP = j - JP
 c      write(*,*) 'JP2=',JP  
1332        a1 = 0.5*(wk(i,jp+1)+wk(i,jp-1)) - wk(i,jp)        a1 = 0.5*(wk(i,jp+1)+wk(i,jp-1)) - wk(i,jp)
1333        b1 = 0.5*(wk(i,jp+1)-wk(i,jp-1))        b1 = 0.5*(wk(i,jp+1)-wk(i,jp-1))
 c      write(*,*) 'a1=',a1,'b1=',b1  
1334        ady(i,j) = wk(i,jp) + rv*(a1*rv + b1) - wk(i,j)        ady(i,j) = wk(i,jp) + rv*(a1*rv + b1) - wk(i,j)
1335        enddo        enddo
1336        enddo        enddo
1337  c      write(*,*) 'toto 2'  
 C  
1338        ELSEIF(IAD.eq.1) then        ELSEIF(IAD.eq.1) then
1339          do j=j1-1,j2+1          do j=j1-1,j2+1
1340        do i=1,imr        do i=1,imr
# Line 1980  C Poles. Line 1961  C Poles.
1961        endif        endif
1962        return        return
1963        end        end
 C  
       subroutine zflip(q,im,km,nc)  
 C This routine flip the array q (in the vertical).  
       real q(im,km,nc)  
 C local dynamic array  
       real qtmp(im,km)  
 C  
       do 4000 IC = 1, nc  
 C  
       do 1000 k=1,km  
       do 1000 i=1,im  
       qtmp(i,k) = q(i,km+1-k,IC)  
 1000  continue  
 C  
       do 2000 i=1,im*km  
 2000  q(i,1,IC) = qtmp(i,1)  
 4000  continue  
       return  
       end  
Legend:
Removed from v.10  
changed lines
  Added in v.32
  ViewVC Help
Powered by ViewVC 1.1.21