compart.f

对工业生产过程结晶过程的一个仿真程序软件包

              ENDDO
           ENDDO 
*          compute f in i-direction            
           DO j=1, meshj3-1
              DO i=1, meshi3-1
                 if((ftemp1(i+1,j)-ftemp1(i,j)).EQ.0.0D0)then
                    i2theta=1.0D10*(ftemp1(i,j)-ftemp1(i-1,j))
                 else
                    i2theta=(ftemp1(i,j)-ftemp1(i-1,j))/
     &                   (ftemp1(i+1,j)-ftemp1(i,j))
                 endif
                 
                 if(((ftemp1(i,j)-ftemp1(i-1,j)).EQ.0.0D0))then
                    i1theta=1.0D10*(ftemp1(i-1,j)-ftemp1(i-2,j))
                 else
                    i1theta=(ftemp1(i-1,j)-ftemp1(i-2,j))/
     &                   (ftemp1(i,j)-ftemp1(i-1,j))
                 endif
                 i1phi=(DABS(i1theta) +i1theta)/(1.0D0+DABS(i1theta))
                 i2phi=(DABS(i2theta) +i2theta)/(1.0D0+DABS(i2theta))
*                High resolution
                 f(i,j)=ftemp1(i,j)-(vei-0.5D0*vei*(1.0D0-vei)
     &                *i1phi)*(ftemp1(i,j)-ftemp1(i-1,j))-0.5D0*vei*
     &                (1.0D0-vei)*i2phi*(ftemp1(i+1,j)-ftemp1(i,j)) 
                 if((f(i,j).lt.-0.1D0).and.(i.lt.5).and.(j.lt.5))then
                    write(myrank+200,150)'f2=', f(i,j),i,j
                 end if       
              ENDDO
           ENDDO      
        ELSE
*          Perfect mixing
*          Compute f in j-direction              
           DO i=1, meshi1
              DO j=1, meshj1              
                 if((f(i,j+1)-f(i,j)).EQ.0.0D0)then
                    j2theta=1.0D10*(f(i,j)-f(i,j-1))
                 else
                    j2theta=(f(i,j)-f(i,j-1))/(f(i,j+1)-f(i,j)) 
                 endif
                 if(((f(i,j)-f(i,j-1)).EQ.0.0D0))then
                    j1theta=1.0D10*(f(i,j-1)-f(i,j-2))
                 else
                    j1theta=(f(i,j-1)-f(i,j-2))/(f(i,j)-f(i,j-1))
                 endif
                 j1phi=(DABS(j1theta) +j1theta)/(1.0D0+DABS(j1theta))
                 j2phi=(DABS(j2theta) +j2theta)/(1.0D0+DABS(j2theta))
                
                 IF (method.EQ.1) THEN
*                   Upwind 
                    ftemp1(i,j)=f(i,j)-vej*(f(i,j)-f(i,j-1)) 
                 ELSE IF (method.EQ.2) THEN
*                   Lax-Wendroff
                    ftemp1(i,j)=f(i,j)-0.5D0*vej*(f(i,j+1)-f(i,j-1))+
     &                   vej*vej*0.5D0*(f(i,j+1)-2.0D0*f(i,j)+f(i,j-1))
                 ELSE
*                   High resolution             
                    ftemp1(i,j)=f(i,j)-(vej-0.5D0*vej*(1.0D0-vej)*j1phi)
     &                   *(f(i,j)-f(i,j-1))-0.5D0*vej*(1.0D0-vej)*j2phi
     &                   *(f(i,j+1)-f(i,j)) 
                 END IF
              ENDDO
           ENDDO
           DO i=meshi2, meshi3-1
              DO j=meshj2, meshj3-1              
                 if((f(i,j+1)-f(i,j)).EQ.0.0D0)then
                    j2theta=1.0D10*(f(i,j)-f(i,j-1))
                 else
                    j2theta=(f(i,j)-f(i,j-1))/(f(i,j+1)-f(i,j)) 
                 endif
                 if(((f(i,j)-f(i,j-1)).EQ.0.0D0))then
                    j1theta=1.0D10*(f(i,j-1)-f(i,j-2))
                 else
                    j1theta=(f(i,j-1)-f(i,j-2))/(f(i,j)-f(i,j-1))
                 endif
                 j1phi=(DABS(j1theta) +j1theta)/(1.0D0+DABS(j1theta))
                 j2phi=(DABS(j2theta) +j2theta)/(1.0D0+DABS(j2theta))         
                 IF (method.EQ.1) THEN
*                   Upwind 
                    ftemp1(i,j)=f(i,j)-vej*(f(i,j)-f(i,j-1)) 
                 ELSE IF (method.EQ.2) THEN
*                   Lax-Wendroff
                    ftemp1(i,j)=f(i,j)-0.5D0*vej*(f(i,j+1)-f(i,j-1))+
     &                   vej*vej*0.5D0*(f(i,j+1)-2.0D0*f(i,j)+f(i,j-1))
                 ELSE
*                   High resolution             
                    ftemp1(i,j)=f(i,j)-(vej-0.5D0*vej*(1.0D0-vej)*j1phi)
     &                   *(f(i,j)-f(i,j-1))-0.5D0*vej*(1.0D0-vej)*j2phi
     &                   *(f(i,j+1)-f(i,j)) 
                 END IF
              ENDDO
           ENDDO
*     compute f in i-direction
           DO j=1, meshj1
              DO i=1, meshi1
                 IF((ftemp1(i+1,j)-ftemp1(i,j)).EQ.0.0D0)THEN
                    i2theta=1.0D10*(ftemp1(i,j)-ftemp1(i-1,j))
                 ELSE
                    i2theta=(ftemp1(i,j)-ftemp1(i-1,j))/
     &                   (ftemp1(i+1,j)-ftemp1(i,j))
                 ENDIF
                 IF(((ftemp1(i,j)-ftemp1(i-1,j)).EQ.0.0D0))THEN
                    i1theta=1.0D10*(ftemp1(i-1,j)-ftemp1(i-2,j))
                 ELSE
                    i1theta=(ftemp1(i-1,j)-ftemp1(i-2,j))/
     &                   (ftemp1(i,j)-ftemp1(i-1,j))
                 ENDIF
                 i1phi=(DABS(i1theta) +i1theta)/(1.0D0+DABS(i1theta))
                 i2phi=(DABS(i2theta) +i2theta)/(1.0D0+DABS(i2theta))
                 IF(method.EQ.1) THEN
*                   Upwind
                    f(i,j)=ftemp1(i,j)-vei*(ftemp1(i,j)-ftemp1(i-1,j))
                 ELSE IF(method.EQ.2) then 
*                   Lax-Wendroff
                    f(i,j)=ftemp1(i,j)-0.5D0*vei*(ftemp1(i+1,j)
     &                   -ftemp1(i-1,j))+vei*vei*0.5D0*(ftemp1(i+1,j)
     &                   -2.0D0*ftemp1(i,j)+ftemp1(i-1,j))
                 ELSE
*                   High resolution
                    f(i,j)=ftemp1(i,j)-(vei-0.5D0*vei*(1.0D0-vei)*i1phi)
     &                *(ftemp1(i,j)-ftemp1(i-1,j))-0.5D0*vei*(1.0D0-vei)
     &                *i2phi*(ftemp1(i+1,j)-ftemp1(i,j))
                 END IF
              ENDDO
           ENDDO
           DO j=meshj2, meshj3-1
              DO i=meshi2, meshi3-1
                 IF((ftemp1(i+1,j)-ftemp1(i,j)).EQ.0.0D0)THEN
                    i2theta=1.0D10*(ftemp1(i,j)-ftemp1(i-1,j))
                 ELSE
                    i2theta=(ftemp1(i,j)-ftemp1(i-1,j))/
     &                   (ftemp1(i+1,j)-ftemp1(i,j))
                 ENDIF
                 IF(((ftemp1(i,j)-ftemp1(i-1,j)).EQ.0.0D0))THEN
                    i1theta=1.0D10*(ftemp1(i-1,j)-ftemp1(i-2,j))
                 ELSE
                    i1theta=(ftemp1(i-1,j)-ftemp1(i-2,j))/
     &                   (ftemp1(i,j)-ftemp1(i-1,j))
                 ENDIF
                 i1phi=(DABS(i1theta) +i1theta)/(1.0D0+DABS(i1theta))
                 i2phi=(DABS(i2theta) +i2theta)/(1.0D0+DABS(i2theta))
                 IF(method.EQ.1) THEN
*                   Upwind
                    f(i,j)=ftemp1(i,j)-vei*(ftemp1(i,j)-ftemp1(i-1,j))
                 ELSE IF(method.EQ.2) then 
*                   Lax-Wendroff
                    f(i,j)=ftemp1(i,j)-0.5D0*vei*(ftemp1(i+1,j)
     &                   -ftemp1(i-1,j))+vei*vei*0.5D0*(ftemp1(i+1,j)
     &                   -2.0D0*ftemp1(i,j)+ftemp1(i-1,j))
                 ELSE
*                   High resolution
                    f(i,j)=ftemp1(i,j)-(vei-0.5D0*vei*(1.0D0-vei)*i1phi)
     &               *(ftemp1(i,j)-ftemp1(i-1,j))-0.5D0*vei*(1.0D0-vei)
     &               *i2phi*(ftemp1(i+1,j)-ftemp1(i,j))
                 END IF
              ENDDO
           ENDDO
        END IF
*       Compute current moments
        IF(numprocs.GT.1)THEN
           DO i=1,meshi3-1
              DO j=1,meshj3-1  
                 mu00(k)=mu00(k)+(f(i,j)+f(i,j+1)+f(i+1,j)+f(i+1,j+1))
     &                /4.0D0*delR1*delR2
                 mu10(k)=mu10(k)
     &                +(f(i,j)*dble(i-1)
     &                +f(i,j+1)*dble(i-1)
     &                +f(i+1,j)*dble(i)
     &                +f(i+1,j+1)*dble(i))/4.0D0
     &                *delR1**2*delR2     
                 mu01(k)=mu01(k)
     &                +(f(i,j)*dble(j-1)
     &                +f(i,j+1)*dble(j)
     &                +f(i+1,j)*dble(j-1)
     &                +f(i+1,j+1)*dble(j))/4.0D0
     &                *delR1*delR2**2
                 mu11(k)=mu11(k)
     &                +(f(i,j)*dble(i-1)*dble(j-1)
     &                +f(i,j+1)*dble(i-1)*dble(j)
     &                +f(i+1,j)*dble(i)*dble(j-1)
     &                +f(i+1,j+1)*dble(i)*dble(j))/4.0D0
     &                *delR1**2*delR2**2     
                 mu20(k)=mu20(k)
     &                +(f(i,j)*dble(i-1)**2
     &                +f(i,j+1)*dble(i-1)**2
     &                +f(i+1,j)*dble(i)**2
     &                +f(i+1,j+1)*dble(i)**2)/4.0D0
     &                *delR1**3*delR2
                 mu02(k)=mu02(k)
     &                +(f(i,j)*dble(j-1)**2
     &                +f(i,j+1)*dble(j)**2
     &                +f(i+1,j)*dble(j-1)**2
     &                +f(i+1,j+1)*dble(j)**2)/4.0D0
     &                *delR1*delR2**3
                 mu21(k)=mu21(k)
     &                +(f(i,j)*dble(i-1)**2*dble(j-1)
     &                +f(i,j+1)*dble(i-1)**2*dble(j)
     &                +f(i+1,j)*dble(i)**2*dble(j-1)
     &                +f(i+1,j+1)*dble(i)**2*dble(j))/4.0D0
     &                *delR1**3*delR2**2
                 mu12(k)=mu12(k)
     &                +(f(i,j)*dble(i-1)*dble(j-1)**2
     &                +f(i,j+1)*dble(i-1)*dble(j)**2
     &                +f(i+1,j)*dble(i)*dble(j-1)**2
     &                +f(i+1,j+1)*dble(i)*dble(j)**2)/4.0D0
     &                *delR1**2*delR2**3
                 mu22(k)=mu22(k)
     &                +(f(i,j)*dble(i-1)**2*dble(j-1)**2
     &                +f(i,j+1)*dble(i-1)**2*dble(j)**2
     &                +f(i+1,j)*dble(i)**2*dble(j-1)**2
     &                +f(i+1,j+1)*dble(i)**2*dble(j)**2)/4.0D0
     &                *delR1**3*delR2**3
                 mu30(k)=mu30(k)
     &                +(f(i,j)*dble(i-1)**3
     &                +f(i,j+1)*dble(i-1)**3
     &                +f(i+1,j)*dble(i)**3
     &                +f(i+1,j+1)*dble(i)**3)/4.0D0
     &                *delR1**4*delR2      
                 mu31(k)=mu31(k)
     &                +(f(i,j)*dble(i-1)**3*dble(j-1)
     &                +f(i,j+1)*dble(i-1)**3*dble(j)
     &                +f(i+1,j)*dble(i)**3*dble(j-1)
     &                +f(i+1,j+1)*dble(i)**3*dble(j))/4.0D0
     &                *delR1**4*delR2**2
              ENDDO
           ENDDO
           volume(k)=mu21(k)-2.0D0/3.0D0*mu30(k)
        ELSE
*          Perfect mixing
           DO i=1,meshi1
              DO j=1,meshj1  
                 mu00(k)=mu00(k)+(f(i,j)+f(i,j+1)+f(i+1,j)+f(i+1,j+1))
     &                /4.0D0*delR1*delR2
                 mu10(k)=mu10(k)
     &                +(f(i,j)*dble(i-1)
     &                +f(i,j+1)*dble(i-1)
     &                +f(i+1,j)*dble(i)
     &                +f(i+1,j+1)*dble(i))/4.0D0
     &                *delR1**2*delR2     
                 mu01(k)=mu01(k)
     &                +(f(i,j)*dble(j-1)
     &                +f(i,j+1)*dble(j)
     &                +f(i+1,j)*dble(j-1)
     &                +f(i+1,j+1)*dble(j))/4.0D0
     &                *delR1*delR2**2
                 mu11(k)=mu11(k)
     &                +(f(i,j)*dble(i-1)*dble(j-1)
     &                +f(i,j+1)*dble(i-1)*dble(j)
     &                +f(i+1,j)*dble(i)*dble(j-1)
     &                +f(i+1,j+1)*dble(i)*dble(j))/4.0D0
     &                *delR1**2*delR2**2     
                 mu20(k)=mu20(k)
     &                +(f(i,j)*dble(i-1)**2
     &                +f(i,j+1)*dble(i-1)**2
     &                +f(i+1,j)*dble(i)**2
     &                +f(i+1,j+1)*dble(i)**2)/4.0D0
     &                *delR1**3*delR2
                 mu02(k)=mu02(k)
     &                +(f(i,j)*dble(j-1)**2
     &                +f(i,j+1)*dble(j)**2
     &                +f(i+1,j)*dble(j-1)**2
     &                +f(i+1,j+1)*dble(j)**2)/4.0D0
     &                *delR1*delR2**3
                 mu21(k)=mu21(k)
     &                +(f(i,j)*dble(i-1)**2*dble(j-1)
     &                +f(i,j+1)*dble(i-1)**2*dble(j)
     &                +f(i+1,j)*dble(i)**2*dble(j-1)
     &                +f(i+1,j+1)*dble(i)**2*dble(j))/4.0D0
     &                *delR1**3*delR2**2
                 mu12(k)=mu12(k)
     &                +(f(i,j)*dble(i-1)*dble(j-1)**2
     &                +f(i,j+1)*dble(i-1)*dble(j)**2
     &                +f(i+1,j)*dble(i)*dble(j-1)**2
     &                +f(i+1,j+1)*dble(i)*dble(j)**2)/4.0D0
     &                *delR1**2*delR2**3
                 mu22(k)=mu22(k)
     &                +(f(i,j)*dble(i-1)**2*dble(j-1)**2
     &                +f(i,j+1)*dble(i-1)**2*dble(j)**2
     &                +f(i+1,j)*dble(i)**2*dble(j-1)**2
     &                +f(i+1,j+1)*dble(i)**2*dble(j)**2)/4.0D0
     &                *delR1**3*delR2**3
                 mu30(k)=mu30(k)
     &                +(f(i,j)*dble(i-1)**3
     &                +f(i,j+1)*dble(i-1)**3
     &                +f(i+1,j)*dble(i)**3
     &                +f(i+1,j+1)*dble(i)**3)/4.0D0
     &                *delR1**4*delR2      
                 mu31(k)=mu31(k)
     &                +(f(i,j)*dble(i-1)**3*dble(j-1)
     &                +f(i,j+1)*dble(i-1)**3*dble(j)
     &                +f(i+1,j)*dble(i)**3*dble(j-1)
     &                +f(i+1,j+1)*dble(i)**3*dble(j))/4.0D0
     &                *delR1**4*delR2**2
              ENDDO
           ENDDO
           DO i=meshi2,meshi3-1
              DO j=meshj2,meshj3-1  
                 mu00(k)=mu00(k)+(f(i,j)+f(i,j+1)+f(i+1,j)+f(i+1,j+1))
     &                /4.0D0*delR1*delR2
                 mu10(k)=mu10(k)
     &                +(f(i,j)*dble(i-1)
     &                +f(i,j+1)*dble(i-1)
     &                +f(i+1,j)*dble(i)
     &                +f(i+1,j+1)*dble(i))/4.0D0
     &                *delR1**2*delR2     
                 mu01(k)=mu01(k)
     &                +(f(i,j)*dble(j-1)
     &                +f(i,j+1)*dble(j)
     &                +f(i+1,j)*dble(j-1)
     &                +f(i+1,j+1)*dble(j))/4.0D0
     &                *delR1*delR2**2
                 mu11(k)=mu11(k)
     &                +(f(i,j)*dble(i-1)*dble(j-1)
     &                +f(i,j+1)*dble(i-1)*dble(j)
     &                +f(i+1,j)*dble(i)*dble(j-1)
     &                +f(i+1,j+1)*dble(i)*dble(j))/4.0D0
     &                *delR1**2*delR2**2     
                 mu20(k)=mu20(k)
     &                +(f(i,j)*dble(i-1)**2
     &                +f(i,j+1)*dble(i-1)**2
     &                +f(i+1,j)*dble(i)**2
     &                +f(i+1,j+1)*dble(i)**2)/4.0D0
     &                *delR1**3*delR2
                 mu02(k)=mu02(k)
     &                +(f(i,j)*dble(j-1)**2
     &                +f(i,j+1)*dble(j)**2
     &                +f(i+1,j)*dble(j-1)**2
     &                +f(i+1,j+1)*dble(j)**2)/4.0D0