demona: D:/demona/2.0/manual/doc/2-31set/problems/test/2/model.f90 Source File

00001 module model
00002 !module modelSet  ! use to differentiate between set and no set models.
00003 
00004 
00005 use global
00006 
00007 use myFunctions
00008 
00009 implicit none
00010 
00011 
00012 contains 
00013 
00014 
00015 function discrete(a)
00016 
00017 
00018 real(8) discrete(globalCurrentSetLength)
00019 
00020 integer i,j,t
00021 
00022 integer a
00023 
00024 integer myI,myJ,myT,myBC,myUnk
00025 
00026 integer myP,myD
00027 
00028 integer,pointer :: IM(:,:,:,:)
00029 
00030 integer,pointer :: IM2(:,:,:,:)
00031 
00032 integer,pointer :: maxI,maxJ,maxT
00033 
00034 integer,pointer :: iVec(:),jVec(:),bcVec(:),tVec(:)
00035 
00036 integer,pointer :: endIndex(:)
00037 
00038 integer modelPhysics
00039 integer modelDomain
00040 
00041 real(8), allocatable :: RHS(:)
00042 
00043 integer,pointer :: dof
00044 
00045 integer,pointer :: nx,ny
00046 
00047 !integer,pointer :: dof
00048 
00049 
00050 ! model specific definitions
00051 
00052 real(8) Tc,Tw,Te,Ts,Tn
00053 
00054 real(8) T1,T2,T3,T4
00055 
00056 real(8) Uc,Uw,Ue,Us,Un
00057 
00058 real(8) U1,U2,U3,U4
00059 
00060 real(8) d2TdX2,d2TdY2
00061 
00062 real(8) d2UdX2,d2UdY2
00063 
00064 
00065 
00066 real(8) dtau
00067 
00068 real(8) upu,upv
00069 
00070 real(8) myU,myV
00071 
00072 
00073 
00074 real(8) sourceT
00075 real(8) sourceU
00076 
00077 real(8) dx,dy
00078 
00079 real(8),pointer :: myLx,myLy
00080 
00081 
00082 
00083 !< if x+eps*v is to be evaluated, then tempy is updated before calling discrete(1)
00084 !< such that tempy(setVec)=tempy(setVec)+lambda*deltax
00085 !< if discrete(0) is called, original y is evaluated
00086 do i=1,a
00087         discreteY=>tempy
00088 end do
00089 
00090 
00091 
00092 
00093 
00094 
00095 !stop
00096 
00097 
00098 
00099 
00100 discrete=-2.d0
00101 
00102 
00103 x=>discreteY
00104 
00105 
00106 ! In this problem, discretization for 1 and 2 is the same although they are located in different
00107 ! part of the domains of different physics
00108 ! they might be coded in a loop so double explicit definition for P1 and P2 is avoided
00109 ! this is not the case for BC 3 since it is defined explicitly
00110 ! however, if and interpolation interface is used, than the bourden is passed to the interface routine
00111 ! so that the discretizations looks the same. After then, P1 and P2 can be defined in a loop.
00112 
00113 
00114 
00115 
00116 
00117 
00118 
00119 
00120 !!!!!!!!!!!!!!
00121 !< assume that  physics-1 starts here
00122 !< also that the domain-1 of physics-1 starts here, as well
00123 !#############################################################################################
00124 !#############################################################################################
00125 !#############################################################################################
00126 !
00127 !############                                   ###
00128 !#############                                  ###
00129 !###            ###                                     ###
00130 !###            ###                                     ###
00131 !###            ###                                     ###
00132 !###            ###                                     ###
00133 !#############          ####            ###
00134 !############           ####            ###
00135 !###                                                    ###
00136 !###                                                    ###
00137 !###                                                    ###
00138 !###                                                    ###
00139 !###                                                    ###
00140 !###                                                    ###     
00141 !
00142 !#############################################################################################
00143 !#############################################################################################
00144 !#############################################################################################
00145 
00146 
00147 !modelPhysics=currentPhysics
00148 !modelDomain=currentDomain
00149 
00150 modelPhysics=1
00151 
00152 sourceT=1.d0
00153 sourceU=-1.d0
00154 
00155 
00156 
00157 !< the idea is  mymodel%discrete%set(1)
00158 !<                              mymodel%discrete%set(2)
00159 !<                              mymodel%discrete%set(3)
00160 ! with set vec pointing to the set of the physics(..)%domain%(..)%set
00161 
00162 
00163 
00164 
00165 
00166 !print*,'dx,dy=',dx,dy
00167 !stop
00168 
00169 
00170 dof=>myModel%physics(modelPhysics)%dof
00171 
00172 
00173 
00174 
00175 allocate(RHS(dof))
00176 RHS=0.d0
00177 
00178 !#####################################################################
00179 !user supplied variables
00180 !#####################################################################
00181 
00182 do myP=1,myModel%endIndex(modelPhysics)
00183 
00184         !#############################################################################################
00185         !
00186         !########                       ##
00187         !##             ##                      ##
00188         !##             ##                      ##
00189         !##             ##      ###             ##
00190         !##             ##                      ##
00191         !##             ##                      ##
00192         !########                       ##
00193         !
00194         !#############################################################################################
00195         
00196         modelDomain=1
00197 
00198         IM=>myModel%physics(modelPhysics)%domain(modelDomain)%dofMatrix
00199 
00200         myLx=>myModel%physics(modelPhysics)%domain(modelDomain)%Lx
00201         myLy=>myModel%physics(modelPhysics)%domain(modelDomain)%Ly
00202 
00203         dx=myLx/(myModel%physics(modelPhysics)%domain(modelDomain)%maxI-myModel%physics(modelPhysics)%domain(modelDomain)%minI)
00204         dy=myLy/(myModel%physics(modelPhysics)%domain(modelDomain)%maxJ-myModel%physics(modelPhysics)%domain(modelDomain)%minJ)
00205         ! normally, especially in finite volume coordinates should be read so that face lengths and center-to-center distances
00206         ! are calculated. for this uniform finite difference grid dx,dy is enough to proceed
00207 
00208         ! caution, how is dx dy in finite volume!!! was lx for the whole domain with ghostcells 
00209         ! or just the real interior domain w/o ghost cells!!! - if so number of ghostcells connection!!
00210 
00211         !print*,dx,dy
00212 
00213 
00214 
00215         do myD=1,myModel%physics(modelPhysics)%endIndex(modelDomain)
00216 
00217 
00218         !boundary
00219 
00220         !---------------------------------------------------------------------------------------------
00221         !all except right
00222         !bc=2
00223         myBC=2   
00224 
00225         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
00226         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
00227 
00228 
00229         do myJ=1,maxJ
00230                 
00231                 j=jVec(myJ)
00232                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
00233                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
00234 
00235                 RHS=0
00236                 do myI=1,maxI
00237 
00238                         i=iVec(myI)
00239                         
00240                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
00241                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
00242                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
00243 
00244                         !unk1
00245                         t=1
00246                         do myT=1,endIndex(t)
00247 
00248                                 ! center values
00249                                 TC=x(IM(i,j,t,absolute)) !< or TC=x(IM(i,j,*1*,absolute))
00250 
00251                                 !print*,TC,i,j
00252 
00253                                 RHS(t)=TC-0.d0
00254 
00255                         end do
00256 
00257                         !unk1
00258                         t=2
00259                         do myT=1,endIndex(t)
00260 
00261                                 ! center values
00262                                 UC=x(IM(i,j,t,absolute))
00263 
00264                                 !print*,TC,i,j
00265 
00266                                 RHS(t)=UC-0.d0
00267 
00268                         end do
00269 
00270 
00271 
00272                         do myT=1,maxT
00273                                 
00274                                 t=tVec(myT)
00275 
00276                                 !print*,IM(i,j,t,relative)
00277 
00278                                 discrete(IM(i,j,t,relative))=rhs(t)
00279 
00280                         end do
00281 
00282                 end do
00283 
00284 
00285 
00286         end do
00287         !---------------------------------------------------------------------------------------------
00288         !stop
00289         !---------------------------------------------------------------------------------------------
00290         
00291         IM2=>myModel%physics(modelPhysics)%domain(2)%dofMatrix
00292         
00293         !right
00294         !bc=3
00295         myBC=3 
00296 
00297         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
00298         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
00299 
00300 
00301         do myJ=1,maxJ
00302                 
00303                 j=jVec(myJ)
00304                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
00305                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
00306 
00307                 RHS=0
00308                 do myI=1,maxI
00309 
00310                         i=iVec(myI)
00311                         
00312                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
00313                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
00314                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
00315 
00316                         !unk1
00317                         t=1
00318                         do myT=1,endIndex(t)
00319 
00320                                 ! center values
00321                                 TC=x(IM(i,j,t,absolute))
00322 
00323                                 !print*,TC,i,j
00324 
00325 
00326                                 ! here the coupling is defined explicitly
00327                                 ! need and interpolation interface
00328                                 ! such as:
00329                                 !
00330                                 ! TC - interpolate(TC_related_info_i,j,t,physics,domain)
00331                         
00332 
00333                                 T1=x(IM2(1,j-1,t,absolute))
00334                                 T2=x(IM2(1+1,j-1,t,absolute))
00335                                 T3=x(IM2(1+1,j,t,absolute))
00336                                 T4=x(IM2(1,j,t,absolute))
00337 
00338                                 !interpolation
00339                                 RHS(t)=TC-(T1+T2+T3+T4)/4.d0
00340                                 !RHS(t)=TC
00341 
00342                         end do
00343 
00344                         !unk2
00345                         t=2
00346                         do myT=1,endIndex(t)
00347 
00348                                 ! center values
00349                                 UC=x(IM(i,j,t,absolute))
00350 
00351                                 !print*,TC,i,j
00352 
00353 
00354                                 ! here the coupling is defined explicitly
00355                                 ! need and interpolation interface
00356                                 ! such as:
00357                                 !
00358                                 ! TC - interpolate(TC_related_info_i,j,t,physics,domain)
00359                         
00360 
00361                                 U1=x(IM2(1,j-1,t,absolute))
00362                                 U2=x(IM2(1+1,j-1,t,absolute))
00363                                 U3=x(IM2(1+1,j,t,absolute))
00364                                 U4=x(IM2(1,j,t,absolute))
00365 
00366                                 !interpolation
00367                                 RHS(t)=UC-(U1+U2+U3+U4)/4.d0
00368                                 !RHS(t)=UC
00369 
00370                         end do
00371 
00372 
00373 
00374 
00375 
00376                         do myT=1,maxT
00377                                 
00378                                 t=tVec(myT)
00379 
00380                                 !print*,IM(i,j,t,relative)
00381 
00382                                 discrete(IM(i,j,t,relative))=rhs(t)
00383 
00384 
00385                         end do
00386 
00387                 end do
00388 
00389 
00390         end do
00391         !---------------------------------------------------------------------------------------------
00392         !stop
00393         !---------------------------------------------------------------------------------------------
00394         !interior
00395         !bc=1
00396         myBC=1
00397 
00398         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
00399         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
00400 
00401         do myJ=1,maxJ
00402                 
00403                 j=jVec(myJ)
00404                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
00405                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
00406 
00407                 RHS=0
00408                 do myI=1,maxI
00409 
00410                         i=iVec(myI)
00411                         
00412                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
00413                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
00414                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
00415 
00416                         !unk1
00417                         t=1
00418                         do myT=1,endIndex(t)
00419 
00420 
00421                                 TC=x(IM(i,j,t,absolute))
00422                                 TW=x(IM(i-1,j,t,absolute))
00423                                 TE=x(IM(i+1,j,t,absolute))
00424                                 TS=x(IM(i,j-1,t,absolute))
00425                                 TN=x(IM(i,j+1,t,absolute))
00426 
00427                                 d2TdX2=(TE+TW-2.d0*TC)/(dx**2.d0)
00428 
00429                                 d2TdY2=(TN+TS-2.d0*TC)/(dy**2.d0)
00430 
00431                                 RHS(t)=d2TdX2+d2TdY2+sourceT
00432 
00433                                 !RHS(t)=TC
00434 
00435                                 !print*,TC,i,j
00436 
00437                         end do
00438 
00439                         !unk1
00440                         t=2
00441                         do myT=1,endIndex(t)
00442 
00443 
00444                                 UC=x(IM(i,j,t,absolute))
00445                                 UW=x(IM(i-1,j,t,absolute))
00446                                 UE=x(IM(i+1,j,t,absolute))
00447                                 US=x(IM(i,j-1,t,absolute))
00448                                 UN=x(IM(i,j+1,t,absolute))
00449 
00450                                 d2UdX2=(UE+UW-2.d0*UC)/(dx**2.d0)
00451 
00452                                 d2UdY2=(UN+US-2.d0*UC)/(dy**2.d0)
00453 
00454                                 RHS(t)=d2UdX2+d2UdY2+sourceU
00455 
00456                                 !RHS(t)=TC
00457 
00458                                 !print*,TC,i,j
00459 
00460                         end do
00461 
00462 
00463 
00464                         do myT=1,maxT
00465 
00466                                 t=tVec(myT)
00467 
00468                                 !print*,IM(i,j,t,relative),rhs(t)
00469 
00470                                 discrete(IM(i,j,t,relative))=rhs(t)
00471 
00472                         end do
00473 
00474                 end do
00475 
00476 
00477         end do
00478         !---------------------------------------------------------------------------------------------
00479         !stop
00480 
00481         end do ! domain 1 loop
00482 
00483 
00484 
00485         !#############################################################################################
00486         !
00487         !########                       ##      ##
00488         !##             ##                      ##      ##
00489         !##             ##                      ##      ##
00490         !##             ##      ###             ##      ##
00491         !##             ##                      ##      ##
00492         !##             ##                      ##      ##
00493         !########                       ##      ##
00494         !
00495         !#############################################################################################
00496 
00497         !start domain 2 loop
00498         modelDomain=2
00499 
00500         IM=>myModel%physics(modelPhysics)%domain(modelDomain)%dofMatrix
00501 
00502         myLx=>myModel%physics(modelPhysics)%domain(modelDomain)%Lx
00503         myLy=>myModel%physics(modelPhysics)%domain(modelDomain)%Ly
00504 
00505         dx=myLx/(myModel%physics(modelPhysics)%domain(modelDomain)%maxI-myModel%physics(modelPhysics)%domain(modelDomain)%minI)
00506         dy=myLy/(myModel%physics(modelPhysics)%domain(modelDomain)%maxJ-myModel%physics(modelPhysics)%domain(modelDomain)%minJ)
00507 
00508         do myD=1,myModel%physics(modelPhysics)%endIndex(modelDomain)
00509 
00510 
00511         !boundary
00512 
00513         !---------------------------------------------------------------------------------------------
00514         !all except left and right
00515         !bc=2
00516         myBC=2   
00517 
00518         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
00519         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
00520 
00521 
00522         do myJ=1,maxJ
00523                 
00524                 j=jVec(myJ)
00525                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
00526                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
00527 
00528                 RHS=0
00529                 do myI=1,maxI
00530 
00531                         i=iVec(myI)
00532                         
00533                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
00534                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
00535                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
00536 
00537                         !unk1
00538                         t=1
00539                         do myT=1,endIndex(t)
00540 
00541                                 ! center values
00542                                 TC=x(IM(i,j,t,absolute))
00543 
00544                                 !print*,TC,i,j
00545 
00546                                 RHS(t)=TC-0.d0
00547 
00548                         end do
00549 
00550                         !unk2
00551                         t=2
00552                         do myT=1,endIndex(t)
00553 
00554                                 ! center values
00555                                 UC=x(IM(i,j,t,absolute))
00556 
00557                                 !print*,TC,i,j
00558 
00559                                 RHS(t)=UC-0.d0
00560 
00561                         end do
00562 
00563 
00564 
00565 
00566                         do myT=1,maxT
00567                                 
00568                                 t=tVec(myT)
00569 
00570                                 !print*,IM(i,j,t,relative)
00571 
00572                                 discrete(IM(i,j,t,relative))=rhs(t)
00573 
00574                         end do
00575 
00576                 end do
00577 
00578 
00579 
00580         end do
00581         !---------------------------------------------------------------------------------------------
00582         !stop
00583         !---------------------------------------------------------------------------------------------
00584 
00585         IM2=>myModel%physics(modelPhysics)%domain(1)%dofMatrix
00586         nx=>myModel%physics(modelPhysics)%domain(1)%nx
00587         ny=>myModel%physics(modelPhysics)%domain(1)%ny
00588 
00589         !left
00590         !bc=3
00591         myBC=3 
00592 
00593         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
00594         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
00595 
00596 
00597         do myJ=1,maxJ
00598                 
00599                 j=jVec(myJ)
00600                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
00601                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
00602 
00603                 RHS=0
00604                 do myI=1,maxI
00605 
00606                         i=iVec(myI)
00607                         
00608                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
00609                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
00610                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
00611 
00612                         !unk1
00613                         t=1
00614                         do myT=1,endIndex(t)
00615 
00616                                 ! center values
00617                                 TC=x(IM(i,j,t,absolute))
00618 
00619                                 !print*,TC,i,j
00620 
00621 
00622                                 ! here the coupling is defined explicitly
00623                                 ! need and interpolation interface
00624                                 ! such as:
00625                                 !
00626                                 ! TC - interpolate(TC_related_info_i,j,t,physics,domain)
00627                         
00628 
00629                                 T1=x(IM2(nx-1,j,t,absolute))
00630                                 T2=x(IM2(nx,j,t,absolute))
00631                                 T3=x(IM2(nx,j+1,t,absolute))
00632                                 T4=x(IM2(nx-1,j+1,t,absolute))
00633 
00634                                 !interpolation
00635                                 RHS(t)=TC-(T1+T2+T3+T4)/4.d0
00636 
00637                         end do
00638 
00639                         !unk1
00640                         t=2
00641                         do myT=1,endIndex(t)
00642 
00643                                 ! center values
00644                                 UC=x(IM(i,j,t,absolute))
00645 
00646                                 !print*,TC,i,j
00647 
00648 
00649                                 ! here the coupling is defined explicitly
00650                                 ! need and interpolation interface
00651                                 ! such as:
00652                                 !
00653                                 ! TC - interpolate(TC_related_info_i,j,t,physics,domain)
00654                         
00655 
00656                                 U1=x(IM2(nx-1,j,t,absolute))
00657                                 U2=x(IM2(nx,j,t,absolute))
00658                                 U3=x(IM2(nx,j+1,t,absolute))
00659                                 U4=x(IM2(nx-1,j+1,t,absolute))
00660 
00661                                 !interpolation
00662                                 RHS(t)=UC-(U1+U2+U3+U4)/4.d0
00663 
00664                         end do
00665 
00666 
00667 
00668                         do myT=1,maxT
00669                                 
00670                                 t=tVec(myT)
00671 
00672                                 !print*,IM(i,j,t,relative)
00673 
00674                                 discrete(IM(i,j,t,relative))=rhs(t)
00675 
00676 
00677                         end do
00678 
00679                 end do
00680 
00681 
00682         end do
00683         !---------------------------------------------------------------------------------------------
00684         !stop
00685         !---------------------------------------------------------------------------------------------
00686 
00687         IM2=>myModel%physics(2)%domain(1)%dofMatrix
00688         
00689         !right
00690         !bc=4
00691         myBC=4 
00692 
00693         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
00694         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
00695 
00696 
00697         do myJ=1,maxJ
00698                 
00699                 j=jVec(myJ)
00700                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
00701                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
00702 
00703                 RHS=0
00704                 do myI=1,maxI
00705 
00706                         i=iVec(myI)
00707                         
00708                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
00709                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
00710                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
00711 
00712                         !unk1
00713                         t=1
00714                         do myT=1,endIndex(t)
00715 
00716                                 ! center values
00717                                 TC=x(IM(i,j,t,absolute))
00718 
00719                                 !print*,TC,i,j
00720 
00721 
00722                                 ! here the coupling is defined explicitly
00723                                 ! need and interpolation interface
00724                                 ! such as:
00725                                 !
00726                                 ! TC - interpolate(TC_related_info_i,j,t,physics,domain)
00727                         
00728 
00729                                 T1=x(IM2(1,j,t,absolute))
00730                                 T2=x(IM2(2,j,t,absolute))
00731                                 T3=x(IM2(2,j+1,t,absolute))
00732                                 T4=x(IM2(1,j+1,t,absolute))
00733 
00734                                 !interpolation
00735                                 RHS(t)=TC-(T1+T2+T3+T4)/4.d0
00736 
00737                         end do
00738 
00739                         !unk1
00740                         t=2
00741                         do myT=1,endIndex(t)
00742 
00743                                 ! center values
00744                                 UC=x(IM(i,j,t,absolute))
00745 
00746                                 !print*,TC,i,j
00747 
00748 
00749                                 ! here the coupling is defined explicitly
00750                                 ! need and interpolation interface
00751                                 ! such as:
00752                                 !
00753                                 ! TC - interpolate(TC_related_info_i,j,t,physics,domain)
00754                         
00755 
00756                                 U1=x(IM2(1,j,t,absolute))
00757                                 U2=x(IM2(2,j,t,absolute))
00758                                 U3=x(IM2(2,j+1,t,absolute))
00759                                 U4=x(IM2(1,j+1,t,absolute))
00760 
00761                                 !interpolation
00762                                 RHS(t)=UC-(U1+U2+U3+U4)/4.d0
00763 
00764                         end do
00765 
00766 
00767 
00768                         do myT=1,maxT
00769                                 
00770                                 t=tVec(myT)
00771 
00772                                 !print*,IM(i,j,t,relative)
00773 
00774                                 discrete(IM(i,j,t,relative))=rhs(t)
00775 
00776 
00777                         end do
00778 
00779                 end do
00780 
00781 
00782         end do
00783         !---------------------------------------------------------------------------------------------  
00784         !stop
00785         !---------------------------------------------------------------------------------------------
00786         ! interior
00787         !bc=1
00788         myBC=1
00789 
00790         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
00791         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
00792 
00793         do myJ=1,maxJ
00794                 
00795                 j=jVec(myJ)
00796                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
00797                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
00798 
00799                 RHS=0
00800                 do myI=1,maxI
00801 
00802                         i=iVec(myI)
00803                         
00804                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
00805                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
00806                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
00807 
00808                         !unk1
00809                         t=1
00810                         do myT=1,endIndex(t)
00811 
00812 
00813                                 TC=x(IM(i,j,t,absolute))
00814                                 TW=x(IM(i-1,j,t,absolute))
00815                                 TE=x(IM(i+1,j,t,absolute))
00816                                 TS=x(IM(i,j-1,t,absolute))
00817                                 TN=x(IM(i,j+1,t,absolute))
00818 
00819                                 d2TdX2=(TE+TW-2.d0*TC)/(dx**2.d0)
00820 
00821                                 d2TdY2=(TN+TS-2.d0*TC)/(dy**2.d0)
00822 
00823                                 RHS(t)=d2TdX2+d2TdY2+sourceT
00824 
00825                                 !RHS(t)=TC
00826 
00827                                 !print*,TC,i,j
00828 
00829                         end do
00830 
00831 
00832                         !unk1
00833                         t=2
00834                         do myT=1,endIndex(t)
00835 
00836 
00837                                 UC=x(IM(i,j,t,absolute))
00838                                 UW=x(IM(i-1,j,t,absolute))
00839                                 UE=x(IM(i+1,j,t,absolute))
00840                                 US=x(IM(i,j-1,t,absolute))
00841                                 UN=x(IM(i,j+1,t,absolute))
00842 
00843                                 d2UdX2=(UE+UW-2.d0*UC)/(dx**2.d0)
00844 
00845                                 d2UdY2=(UN+US-2.d0*UC)/(dy**2.d0)
00846 
00847                                 RHS(t)=d2UdX2+d2UdY2+sourceU
00848 
00849                                 !RHS(t)=TC
00850 
00851                                 !print*,TC,i,j
00852 
00853                         end do
00854 
00855 
00856 
00857                         do myT=1,maxT
00858 
00859                                 t=tVec(myT)
00860 
00861                                 !print*,IM(i,j,t,relative),rhs(t)
00862 
00863                                 discrete(IM(i,j,t,relative))=rhs(t)
00864 
00865                         end do
00866 
00867                 end do
00868 
00869 
00870         end do
00871         !---------------------------------------------------------------------------------------------
00872         !stop
00873 
00874         end do ! domain 2 loop
00875 
00876 end do ! physics loop
00877 
00878 !#####################################################################
00879 
00880 
00881 deallocate(RHS)
00882 
00883 !x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x
00884 !x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x
00885 !x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x
00886 !x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x
00887 
00888 !#############################################################################################
00889 !#############################################################################################
00890 !#############################################################################################
00891 !
00892 !############                                   ###             ###
00893 !#############                                  ###             ###
00894 !###            ###                                     ###             ###
00895 !###            ###                                     ###             ###
00896 !###            ###                                     ###             ###
00897 !###            ###                                     ###             ###
00898 !#############          ####            ###             ###
00899 !############           ####            ###             ###
00900 !###                                                    ###             ###
00901 !###                                                    ###             ###
00902 !###                                                    ###             ###
00903 !###                                                    ###             ###
00904 !###                                                    ###             ###
00905 !###                                                    ###             ###
00906 !
00907 !#############################################################################################
00908 !#############################################################################################
00909 !#############################################################################################
00910 
00911 
00912 modelPhysics=2
00913 
00914 
00915 dof=>myModel%physics(modelPhysics)%dof
00916 
00917 
00918 
00919 allocate(RHS(dof))
00920 RHS=0.d0
00921 
00922 !#####################################################################
00923 !user supplied variables
00924 !#####################################################################
00925 
00926 do myP=1,myModel%endIndex(modelPhysics)
00927 
00928         !#############################################################################################
00929         !
00930         !########                       ##
00931         !##             ##                      ##
00932         !##             ##                      ##
00933         !##             ##      ###             ##
00934         !##             ##                      ##
00935         !##             ##                      ##
00936         !########                       ##
00937         !
00938         !#############################################################################################
00939         
00940         modelDomain=1
00941 
00942         IM=>myModel%physics(modelPhysics)%domain(modelDomain)%dofMatrix
00943 
00944         myLx=>myModel%physics(modelPhysics)%domain(modelDomain)%Lx
00945         myLy=>myModel%physics(modelPhysics)%domain(modelDomain)%Ly
00946 
00947         dx=myLx/(myModel%physics(modelPhysics)%domain(modelDomain)%maxI-myModel%physics(modelPhysics)%domain(modelDomain)%minI)
00948         dy=myLy/(myModel%physics(modelPhysics)%domain(modelDomain)%maxJ-myModel%physics(modelPhysics)%domain(modelDomain)%minJ)
00949         ! normally, especially in finite volume coordinates should be read so that face lengths and center-to-center distances
00950         ! are calculated. for this uniform finite difference grid dx,dy is enough to proceed
00951 
00952         ! caution, how is dx dy in finite volume!!! was lx for the whole domain with ghostcells 
00953         ! or just the real interior domain w/o ghost cells!!! - if so number of ghostcells connection!!
00954 
00955 
00956         do myD=1,myModel%physics(modelPhysics)%endIndex(modelDomain)
00957 
00958 
00959         !boundary
00960 
00961         !---------------------------------------------------------------------------------------------
00962         !all except left and right
00963         !bc=2
00964         myBC=2   
00965 
00966         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
00967         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
00968 
00969 
00970         do myJ=1,maxJ
00971                 
00972                 j=jVec(myJ)
00973                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
00974                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
00975 
00976                 RHS=0
00977                 do myI=1,maxI
00978 
00979                         i=iVec(myI)
00980                         
00981                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
00982                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
00983                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
00984 
00985                         !unk1
00986                         t=1
00987                         do myT=1,endIndex(t)
00988 
00989                                 ! center values
00990                                 TC=x(IM(i,j,t,absolute)) !< or TC=x(IM(i,j,*1*,absolute))
00991 
00992                                 !print*,TC,i,j
00993 
00994                                 RHS(t)=TC-0.d0
00995 
00996                         end do
00997 
00998                         !unk1
00999                         t=2
01000                         do myT=1,endIndex(t)
01001 
01002                                 ! center values
01003                                 UC=x(IM(i,j,t,absolute))
01004 
01005                                 !print*,TC,i,j
01006 
01007                                 RHS(t)=UC-0.d0
01008 
01009                         end do
01010 
01011 
01012 
01013                         do myT=1,maxT
01014                                 
01015                                 t=tVec(myT)
01016 
01017                                 !print*,IM(i,j,t,relative)
01018 
01019                                 discrete(IM(i,j,t,relative))=rhs(t)
01020 
01021                         end do
01022 
01023                 end do
01024 
01025 
01026 
01027         end do
01028         !---------------------------------------------------------------------------------------------
01029         !stop
01030         !---------------------------------------------------------------------------------------------
01031 
01032         IM2=>myModel%physics(1)%domain(2)%dofMatrix
01033         nx=>myModel%physics(1)%domain(2)%nx
01034         ny=>myModel%physics(1)%domain(2)%ny
01035 
01036         !left
01037         !bc=3
01038         myBC=3 
01039 
01040         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
01041         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
01042 
01043 
01044         do myJ=1,maxJ
01045                 
01046                 j=jVec(myJ)
01047                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
01048                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
01049 
01050                 RHS=0
01051                 do myI=1,maxI
01052 
01053                         i=iVec(myI)
01054                         
01055                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
01056                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
01057                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
01058 
01059                         !unk1
01060                         t=1
01061                         do myT=1,endIndex(t)
01062 
01063                                 ! center values
01064                                 TC=x(IM(i,j,t,absolute))
01065 
01066                                 !print*,TC,i,j
01067 
01068 
01069                                 ! here the coupling is defined explicitly
01070                                 ! need and interpolation interface
01071                                 ! such as:
01072                                 !
01073                                 ! TC - interpolate(TC_related_info_i,j,t,physics,domain)
01074                         
01075 
01076                                 T1=x(IM2(nx-1,j-1,t,absolute))
01077                                 T2=x(IM2(nx,j-1,t,absolute))
01078                                 T3=x(IM2(nx,j,t,absolute))
01079                                 T4=x(IM2(nx-1,j,t,absolute))
01080 
01081                                 !interpolation
01082                                 RHS(t)=TC-(T1+T2+T3+T4)/4.d0
01083 
01084                         end do
01085 
01086                         !unk2
01087                         t=2
01088                         do myT=1,endIndex(t)
01089 
01090                                 ! center values
01091                                 UC=x(IM(i,j,t,absolute))
01092 
01093                                 !print*,TC,i,j
01094 
01095 
01096                                 ! here the coupling is defined explicitly
01097                                 ! need and interpolation interface
01098                                 ! such as:
01099                                 !
01100                                 ! TC - interpolate(TC_related_info_i,j,t,physics,domain)
01101                         
01102 
01103                                 U1=x(IM2(nx-1,j-1,t,absolute))
01104                                 U2=x(IM2(nx,j-1,t,absolute))
01105                                 U3=x(IM2(nx,j,t,absolute))
01106                                 U4=x(IM2(nx-1,j,t,absolute))
01107 
01108                                 !interpolation
01109                                 RHS(t)=UC-(U1+U2+U3+U4)/4.d0
01110 
01111                         end do
01112 
01113 
01114                         do myT=1,maxT
01115                                 
01116                                 t=tVec(myT)
01117 
01118                                 !print*,IM(i,j,t,relative)
01119 
01120                                 discrete(IM(i,j,t,relative))=rhs(t)
01121 
01122 
01123                         end do
01124 
01125                 end do
01126 
01127 
01128         end do
01129         !---------------------------------------------------------------------------------------------
01130         !stop
01131         !---------------------------------------------------------------------------------------------
01132 
01133         IM2=>myModel%physics(3)%domain(1)%dofMatrix
01134 
01135         !right
01136         !bc=4
01137         myBC=4
01138 
01139         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
01140         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
01141 
01142 
01143         do myJ=1,maxJ
01144                 
01145                 j=jVec(myJ)
01146                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
01147                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
01148 
01149                 RHS=0
01150                 do myI=1,maxI
01151 
01152                         i=iVec(myI)
01153                         
01154                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
01155                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
01156                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
01157 
01158                         !unk1
01159                         t=1
01160                         do myT=1,endIndex(t)
01161 
01162                                 ! center values
01163                                 TC=x(IM(i,j,t,absolute))
01164 
01165                                 !print*,TC,i,j
01166 
01167 
01168                                 ! here the coupling is defined explicitly
01169                                 ! need and interpolation interface
01170                                 ! such as:
01171                                 !
01172                                 ! TC - interpolate(TC_related_info_i,j,t,physics,domain)
01173                         
01174 
01175                                 T1=x(IM2(1,j-1,t,absolute))
01176                                 T2=x(IM2(2,j-1,t,absolute))
01177                                 T3=x(IM2(2,j,t,absolute))
01178                                 T4=x(IM2(1,j,t,absolute))
01179 
01180                                 !interpolation
01181                                 RHS(t)=TC-(T1+T2+T3+T4)/4.d0
01182 
01183                         end do
01184 
01185                         !unk2
01186                         t=2
01187                         do myT=1,endIndex(t)
01188 
01189                                 ! center values
01190                                 UC=x(IM(i,j,t,absolute))
01191 
01192                                 !print*,TC,i,j
01193 
01194 
01195                                 ! here the coupling is defined explicitly
01196                                 ! need and interpolation interface
01197                                 ! such as:
01198                                 !
01199                                 ! TC - interpolate(TC_related_info_i,j,t,physics,domain)
01200                         
01201 
01202                                 U1=x(IM2(1,j-1,t,absolute))
01203                                 U2=x(IM2(2,j-1,t,absolute))
01204                                 U3=x(IM2(2,j,t,absolute))
01205                                 U4=x(IM2(1,j,t,absolute))
01206 
01207                                 !interpolation
01208                                 RHS(t)=UC-(U1+U2+U3+U4)/4.d0
01209 
01210                         end do
01211 
01212 
01213                         do myT=1,maxT
01214                                 
01215                                 t=tVec(myT)
01216 
01217                                 !print*,IM(i,j,t,relative)
01218 
01219                                 discrete(IM(i,j,t,relative))=rhs(t)
01220 
01221 
01222                         end do
01223 
01224                 end do
01225 
01226 
01227         end do
01228         !---------------------------------------------------------------------------------------------
01229         !stop
01230         !---------------------------------------------------------------------------------------------
01231         !interior
01232         !bc=1
01233         myBC=1
01234 
01235         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
01236         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
01237 
01238         do myJ=1,maxJ
01239                 
01240                 j=jVec(myJ)
01241                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
01242                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
01243 
01244                 RHS=0
01245                 do myI=1,maxI
01246 
01247                         i=iVec(myI)
01248                         
01249                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
01250                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
01251                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
01252 
01253                         !unk1
01254                         t=1
01255                         do myT=1,endIndex(t)
01256 
01257 
01258                                 TC=x(IM(i,j,t,absolute))
01259                                 TW=x(IM(i-1,j,t,absolute))
01260                                 TE=x(IM(i+1,j,t,absolute))
01261                                 TS=x(IM(i,j-1,t,absolute))
01262                                 TN=x(IM(i,j+1,t,absolute))
01263 
01264                                 d2TdX2=(TE+TW-2.d0*TC)/(dx**2.d0)
01265 
01266                                 d2TdY2=(TN+TS-2.d0*TC)/(dy**2.d0)
01267 
01268                                 RHS(t)=d2TdX2+d2TdY2+sourceT
01269 
01270                                 !RHS(t)=TC
01271 
01272                                 !print*,TC,i,j
01273 
01274                         end do
01275 
01276                         !unk1
01277                         t=2
01278                         do myT=1,endIndex(t)
01279 
01280 
01281                                 UC=x(IM(i,j,t,absolute))
01282                                 UW=x(IM(i-1,j,t,absolute))
01283                                 UE=x(IM(i+1,j,t,absolute))
01284                                 US=x(IM(i,j-1,t,absolute))
01285                                 UN=x(IM(i,j+1,t,absolute))
01286 
01287                                 d2UdX2=(UE+UW-2.d0*UC)/(dx**2.d0)
01288 
01289                                 d2UdY2=(UN+US-2.d0*UC)/(dy**2.d0)
01290 
01291                                 RHS(t)=d2UdX2+d2UdY2+sourceU
01292 
01293                                 !RHS(t)=TC
01294 
01295                                 !print*,TC,i,j
01296 
01297                         end do
01298 
01299 
01300 
01301                         do myT=1,maxT
01302 
01303                                 t=tVec(myT)
01304 
01305                                 !print*,IM(i,j,t,relative),rhs(t)
01306 
01307                                 discrete(IM(i,j,t,relative))=rhs(t)
01308 
01309                         end do
01310 
01311                 end do
01312 
01313 
01314         end do
01315         !---------------------------------------------------------------------------------------------
01316         !stop
01317 
01318         end do ! domain 1 loop
01319 
01320 end do ! physics loop
01321 
01322 
01323 !#####################################################################
01324 
01325 
01326 
01327 
01328 deallocate(RHS)
01329 
01330 !x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x
01331 !x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x
01332 !x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x
01333 !x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x
01334 
01335 
01336 
01337 
01338 
01339 
01340 !#############################################################################################
01341 !#############################################################################################
01342 !#############################################################################################
01343 !
01344 !############                                   ###             ###             ###
01345 !#############                                  ###             ###             ###
01346 !###            ###                                     ###             ###             ###
01347 !###            ###                                     ###             ###             ###
01348 !###            ###                                     ###             ###             ###
01349 !###            ###                                     ###             ###             ###
01350 !#############          ####            ###             ###             ###
01351 !############           ####            ###             ###             ###
01352 !###                                                    ###             ###             ###
01353 !###                                                    ###             ###             ###
01354 !###                                                    ###             ###             ###
01355 !###                                                    ###             ###             ###
01356 !###                                                    ###             ###             ###
01357 !###                                                    ###             ###             ###
01358 !
01359 !#############################################################################################
01360 !#############################################################################################
01361 !#############################################################################################
01362 
01363 
01364 !modelPhysics=currentPhysics
01365 !modelDomain=currentDomain
01366 
01367 modelPhysics=3
01368 sourceT=1.d0
01369 sourceU=-1.d0
01370 
01371 
01372 
01373 !< the idea is  mymodel%discrete%set(1)
01374 !<                              mymodel%discrete%set(2)
01375 !<                              mymodel%discrete%set(3)
01376 ! with set vec pointing to the set of the physics(..)%domain%(..)%set
01377 
01378 
01379 
01380 
01381 
01382 !print*,'dx,dy=',dx,dy
01383 !stop
01384 
01385 
01386 dof=>myModel%physics(modelPhysics)%dof
01387 
01388 
01389 allocate(RHS(dof))
01390 RHS=0.d0
01391 
01392 !#####################################################################
01393 !user supplied variables
01394 !#####################################################################
01395 
01396 do myP=1,myModel%endIndex(modelPhysics)
01397 
01398         !#############################################################################################
01399         !
01400         !########                       ##
01401         !##             ##                      ##
01402         !##             ##                      ##
01403         !##             ##      ###             ##
01404         !##             ##                      ##
01405         !##             ##                      ##
01406         !########                       ##
01407         !
01408         !#############################################################################################
01409         
01410         !start domain 1 loop
01411         modelDomain=1
01412 
01413         IM=>myModel%physics(modelPhysics)%domain(modelDomain)%dofMatrix
01414 
01415         IM2=>myModel%physics(modelPhysics)%domain(2)%dofMatrix
01416 
01417         myLx=>myModel%physics(modelPhysics)%domain(modelDomain)%Lx
01418         myLy=>myModel%physics(modelPhysics)%domain(modelDomain)%Ly
01419 
01420         dx=myLx/(myModel%physics(modelPhysics)%domain(modelDomain)%maxI-myModel%physics(modelPhysics)%domain(modelDomain)%minI)
01421         dy=myLy/(myModel%physics(modelPhysics)%domain(modelDomain)%maxJ-myModel%physics(modelPhysics)%domain(modelDomain)%minJ)
01422 
01423         do myD=1,myModel%physics(modelPhysics)%endIndex(modelDomain)
01424 
01425 
01426         !boundary
01427 
01428         !---------------------------------------------------------------------------------------------
01429         !all except left and right
01430         !bc=2
01431         myBC=2   
01432 
01433         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
01434         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
01435 
01436 
01437         do myJ=1,maxJ
01438                 
01439                 j=jVec(myJ)
01440                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
01441                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
01442 
01443                 RHS=0
01444                 do myI=1,maxI
01445 
01446                         i=iVec(myI)
01447                         
01448                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
01449                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
01450                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
01451 
01452                         !unk1
01453                         t=1
01454                         do myT=1,endIndex(t)
01455 
01456                                 ! center values
01457                                 TC=x(IM(i,j,t,absolute))
01458 
01459                                 !print*,TC,i,j
01460 
01461                                 RHS(t)=TC-0.d0
01462 
01463                         end do
01464 
01465                         !unk2
01466                         t=2
01467                         do myT=1,endIndex(t)
01468 
01469                                 ! center values
01470                                 UC=x(IM(i,j,t,absolute))
01471 
01472                                 !print*,TC,i,j
01473 
01474                                 RHS(t)=UC-0.d0
01475 
01476                         end do
01477 
01478 
01479 
01480 
01481                         do myT=1,maxT
01482                                 
01483                                 t=tVec(myT)
01484 
01485                                 !print*,IM(i,j,t,relative)
01486 
01487                                 discrete(IM(i,j,t,relative))=rhs(t)
01488 
01489                         end do
01490 
01491                 end do
01492 
01493 
01494 
01495         end do
01496         !---------------------------------------------------------------------------------------------
01497         !stop
01498         !---------------------------------------------------------------------------------------------
01499 
01500         IM2=>myModel%physics(2)%domain(1)%dofMatrix
01501         nx=>myModel%physics(2)%domain(1)%nx
01502         ny=>myModel%physics(2)%domain(1)%ny
01503 
01504         !left
01505         !bc=3
01506         myBC=3 
01507 
01508         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
01509         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
01510 
01511 
01512         do myJ=1,maxJ
01513                 
01514                 j=jVec(myJ)
01515                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
01516                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
01517 
01518                 RHS=0
01519                 do myI=1,maxI
01520 
01521                         i=iVec(myI)
01522                         
01523                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
01524                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
01525                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
01526 
01527                         !unk1
01528                         t=1
01529                         do myT=1,endIndex(t)
01530 
01531                                 ! center values
01532                                 TC=x(IM(i,j,t,absolute))
01533 
01534                                 !print*,TC,i,j
01535 
01536 
01537                                 ! here the coupling is defined explicitly
01538                                 ! need and interpolation interface
01539                                 ! such as:
01540                                 !
01541                                 ! TC - interpolate(TC_related_info_i,j,t,physics,domain)
01542                         
01543 
01544                                 T1=x(IM2(nx-1,j,t,absolute))
01545                                 T2=x(IM2(nx,j,t,absolute))
01546                                 T3=x(IM2(nx,j+1,t,absolute))
01547                                 T4=x(IM2(nx-1,j+1,t,absolute))
01548 
01549                                 !interpolation
01550                                 RHS(t)=TC-(T1+T2+T3+T4)/4.d0
01551 
01552                         end do
01553 
01554                         !unk1
01555                         t=2
01556                         do myT=1,endIndex(t)
01557 
01558                                 ! center values
01559                                 UC=x(IM(i,j,t,absolute))
01560 
01561                                 !print*,TC,i,j
01562 
01563 
01564                                 ! here the coupling is defined explicitly
01565                                 ! need and interpolation interface
01566                                 ! such as:
01567                                 !
01568                                 ! TC - interpolate(TC_related_info_i,j,t,physics,domain)
01569                         
01570 
01571                                 U1=x(IM2(nx-1,j,t,absolute))
01572                                 U2=x(IM2(nx,j,t,absolute))
01573                                 U3=x(IM2(nx,j+1,t,absolute))
01574                                 U4=x(IM2(nx-1,j+1,t,absolute))
01575 
01576                                 !interpolation
01577                                 RHS(t)=UC-(U1+U2+U3+U4)/4.d0
01578 
01579                         end do
01580 
01581 
01582 
01583                         do myT=1,maxT
01584                                 
01585                                 t=tVec(myT)
01586 
01587                                 !print*,IM(i,j,t,relative)
01588 
01589                                 discrete(IM(i,j,t,relative))=rhs(t)
01590 
01591 
01592                         end do
01593 
01594                 end do
01595 
01596 
01597         end do
01598         !---------------------------------------------------------------------------------------------
01599         !stop
01600         !---------------------------------------------------------------------------------------------
01601 
01602         IM2=>myModel%physics(modelPhysics)%domain(2)%dofMatrix
01603         
01604         !right
01605         !bc=4
01606         myBC=4 
01607 
01608         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
01609         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
01610 
01611 
01612         do myJ=1,maxJ
01613                 
01614                 j=jVec(myJ)
01615                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
01616                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
01617 
01618                 RHS=0
01619                 do myI=1,maxI
01620 
01621                         i=iVec(myI)
01622                         
01623                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
01624                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
01625                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
01626 
01627                         !unk1
01628                         t=1
01629                         do myT=1,endIndex(t)
01630 
01631                                 ! center values
01632                                 TC=x(IM(i,j,t,absolute))
01633 
01634                                 !print*,TC,i,j
01635 
01636 
01637                                 ! here the coupling is defined explicitly
01638                                 ! need and interpolation interface
01639                                 ! such as:
01640                                 !
01641                                 ! TC - interpolate(TC_related_info_i,j,t,physics,domain)
01642                         
01643 
01644                                 T1=x(IM2(1,j,t,absolute))
01645                                 T2=x(IM2(2,j,t,absolute))
01646                                 T3=x(IM2(2,j+1,t,absolute))
01647                                 T4=x(IM2(1,j+1,t,absolute))
01648 
01649                                 !interpolation
01650                                 RHS(t)=TC-(T1+T2+T3+T4)/4.d0
01651 
01652                         end do
01653 
01654                         !unk1
01655                         t=2
01656                         do myT=1,endIndex(t)
01657 
01658                                 ! center values
01659                                 UC=x(IM(i,j,t,absolute))
01660 
01661                                 !print*,TC,i,j
01662 
01663 
01664                                 ! here the coupling is defined explicitly
01665                                 ! need and interpolation interface
01666                                 ! such as:
01667                                 !
01668                                 ! TC - interpolate(TC_related_info_i,j,t,physics,domain)
01669                         
01670 
01671                                 U1=x(IM2(1,j,t,absolute))
01672                                 U2=x(IM2(2,j,t,absolute))
01673                                 U3=x(IM2(2,j+1,t,absolute))
01674                                 U4=x(IM2(1,j+1,t,absolute))
01675 
01676                                 !interpolation
01677                                 RHS(t)=UC-(U1+U2+U3+U4)/4.d0
01678 
01679                         end do
01680 
01681 
01682 
01683                         do myT=1,maxT
01684                                 
01685                                 t=tVec(myT)
01686 
01687                                 !print*,IM(i,j,t,relative)
01688 
01689                                 discrete(IM(i,j,t,relative))=rhs(t)
01690 
01691 
01692                         end do
01693 
01694                 end do
01695 
01696 
01697         end do
01698         !---------------------------------------------------------------------------------------------  
01699         !stop
01700         !---------------------------------------------------------------------------------------------
01701         ! interior
01702         !bc=1
01703         myBC=1
01704 
01705         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
01706         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
01707 
01708         do myJ=1,maxJ
01709                 
01710                 j=jVec(myJ)
01711                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
01712                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
01713 
01714                 RHS=0
01715                 do myI=1,maxI
01716 
01717                         i=iVec(myI)
01718                         
01719                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
01720                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
01721                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
01722 
01723                         !unk1
01724                         t=1
01725                         do myT=1,endIndex(t)
01726 
01727 
01728                                 TC=x(IM(i,j,t,absolute))
01729                                 TW=x(IM(i-1,j,t,absolute))
01730                                 TE=x(IM(i+1,j,t,absolute))
01731                                 TS=x(IM(i,j-1,t,absolute))
01732                                 TN=x(IM(i,j+1,t,absolute))
01733 
01734                                 d2TdX2=(TE+TW-2.d0*TC)/(dx**2.d0)
01735 
01736                                 d2TdY2=(TN+TS-2.d0*TC)/(dy**2.d0)
01737 
01738                                 RHS(t)=d2TdX2+d2TdY2+sourceT
01739 
01740                                 !RHS(t)=TC
01741 
01742                                 !print*,TC,i,j
01743 
01744                         end do
01745 
01746 
01747                         !unk1
01748                         t=2
01749                         do myT=1,endIndex(t)
01750 
01751 
01752                                 UC=x(IM(i,j,t,absolute))
01753                                 UW=x(IM(i-1,j,t,absolute))
01754                                 UE=x(IM(i+1,j,t,absolute))
01755                                 US=x(IM(i,j-1,t,absolute))
01756                                 UN=x(IM(i,j+1,t,absolute))
01757 
01758                                 d2UdX2=(UE+UW-2.d0*UC)/(dx**2.d0)
01759 
01760                                 d2UdY2=(UN+US-2.d0*UC)/(dy**2.d0)
01761 
01762                                 RHS(t)=d2UdX2+d2UdY2+sourceU
01763 
01764                                 !RHS(t)=TC
01765 
01766                                 !print*,TC,i,j
01767 
01768                         end do
01769 
01770 
01771 
01772                         do myT=1,maxT
01773 
01774                                 t=tVec(myT)
01775 
01776                                 !print*,IM(i,j,t,relative),rhs(t)
01777 
01778                                 discrete(IM(i,j,t,relative))=rhs(t)
01779 
01780                         end do
01781 
01782                 end do
01783 
01784 
01785         end do
01786         !---------------------------------------------------------------------------------------------
01787         !stop
01788 
01789         end do ! domain 1 loop
01790 
01791 
01792         !#############################################################################################
01793         !
01794         !########                       ##      ##
01795         !##             ##                      ##      ##
01796         !##             ##                      ##      ##
01797         !##             ##      ###             ##      ##
01798         !##             ##                      ##      ##
01799         !##             ##                      ##      ##
01800         !########                       ##      ##
01801         !
01802         !#############################################################################################
01803 
01804 
01805         modelDomain=2
01806 
01807         IM=>myModel%physics(modelPhysics)%domain(modelDomain)%dofMatrix
01808 
01809         myLx=>myModel%physics(modelPhysics)%domain(modelDomain)%Lx
01810         myLy=>myModel%physics(modelPhysics)%domain(modelDomain)%Ly
01811 
01812         dx=myLx/(myModel%physics(modelPhysics)%domain(modelDomain)%maxI-myModel%physics(modelPhysics)%domain(modelDomain)%minI)
01813         dy=myLy/(myModel%physics(modelPhysics)%domain(modelDomain)%maxJ-myModel%physics(modelPhysics)%domain(modelDomain)%minJ)
01814         ! normally, especially in finite volume coordinates should be read so that face lengths and center-to-center distances
01815         ! are calculated. for this uniform finite difference grid dx,dy is enough to proceed
01816 
01817         ! caution, how is dx dy in finite volume!!! was lx for the whole domain with ghostcells 
01818         ! or just the real interior domain w/o ghost cells!!! - if so number of ghostcells connection!!
01819 
01820 
01821         do myD=1,myModel%physics(modelPhysics)%endIndex(modelDomain)
01822 
01823 
01824         !boundary
01825 
01826         !---------------------------------------------------------------------------------------------
01827         !all except left
01828         !bc=2
01829         myBC=2   
01830 
01831         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
01832         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
01833 
01834 
01835         do myJ=1,maxJ
01836                 
01837                 j=jVec(myJ)
01838                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
01839                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
01840 
01841                 RHS=0
01842                 do myI=1,maxI
01843 
01844                         i=iVec(myI)
01845                         
01846                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
01847                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
01848                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
01849 
01850                         !unk1
01851                         t=1
01852                         do myT=1,endIndex(t)
01853 
01854                                 ! center values
01855                                 TC=x(IM(i,j,t,absolute)) !< or TC=x(IM(i,j,*1*,absolute))
01856 
01857                                 !print*,TC,i,j
01858 
01859                                 RHS(t)=TC-0.d0
01860 
01861                         end do
01862 
01863                         !unk1
01864                         t=2
01865                         do myT=1,endIndex(t)
01866 
01867                                 ! center values
01868                                 UC=x(IM(i,j,t,absolute))
01869 
01870                                 !print*,TC,i,j
01871 
01872                                 RHS(t)=UC-0.d0
01873 
01874                         end do
01875 
01876 
01877 
01878                         do myT=1,maxT
01879                                 
01880                                 t=tVec(myT)
01881 
01882                                 !print*,IM(i,j,t,relative)
01883 
01884                                 discrete(IM(i,j,t,relative))=rhs(t)
01885 
01886                         end do
01887 
01888                 end do
01889 
01890 
01891 
01892         end do
01893         !---------------------------------------------------------------------------------------------
01894         !stop
01895         !---------------------------------------------------------------------------------------------
01896         !right
01897         !bc=3
01898         myBC=3 
01899 
01900         IM2=>myModel%physics(modelPhysics)%domain(1)%dofMatrix
01901         nx=>myModel%physics(modelPhysics)%domain(1)%nx
01902         ny=>myModel%physics(modelPhysics)%domain(1)%ny
01903 
01904         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
01905         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
01906 
01907 
01908         do myJ=1,maxJ
01909                 
01910                 j=jVec(myJ)
01911                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
01912                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
01913 
01914                 RHS=0
01915                 do myI=1,maxI
01916 
01917                         i=iVec(myI)
01918                         
01919                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
01920                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
01921                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
01922 
01923                         !unk1
01924                         t=1
01925                         do myT=1,endIndex(t)
01926 
01927                                 ! center values
01928                                 TC=x(IM(i,j,t,absolute))
01929 
01930                                 !print*,TC,i,j
01931 
01932 
01933                                 ! here the coupling is defined explicitly
01934                                 ! need and interpolation interface
01935                                 ! such as:
01936                                 !
01937                                 ! TC - interpolate(TC_related_info_i,j,t,physics,domain)
01938                         
01939 
01940                                 T1=x(IM2(nx-1,j-1,t,absolute))
01941                                 T2=x(IM2(nx,j-1,t,absolute))
01942                                 T3=x(IM2(nx,j,t,absolute))
01943                                 T4=x(IM2(nx-1,j,t,absolute))
01944 
01945                                 !interpolation
01946                                 RHS(t)=TC-(T1+T2+T3+T4)/4.d0
01947 
01948                         end do
01949 
01950                         !unk2
01951                         t=2
01952                         do myT=1,endIndex(t)
01953 
01954                                 ! center values
01955                                 UC=x(IM(i,j,t,absolute))
01956 
01957                                 !print*,TC,i,j
01958 
01959 
01960                                 ! here the coupling is defined explicitly
01961                                 ! need and interpolation interface
01962                                 ! such as:
01963                                 !
01964                                 ! TC - interpolate(TC_related_info_i,j,t,physics,domain)
01965                         
01966 
01967                                 U1=x(IM2(nx-1,j-1,t,absolute))
01968                                 U2=x(IM2(nx,j-1,t,absolute))
01969                                 U3=x(IM2(nx,j,t,absolute))
01970                                 U4=x(IM2(nx-1,j,t,absolute))
01971 
01972                                 !interpolation
01973                                 RHS(t)=UC-(U1+U2+U3+U4)/4.d0
01974 
01975                         end do
01976 
01977 
01978 
01979 
01980 
01981                         do myT=1,maxT
01982                                 
01983                                 t=tVec(myT)
01984 
01985                                 !print*,IM(i,j,t,relative)
01986 
01987                                 discrete(IM(i,j,t,relative))=rhs(t)
01988 
01989 
01990                         end do
01991 
01992                 end do
01993 
01994 
01995         end do
01996         !---------------------------------------------------------------------------------------------
01997         !stop
01998         !---------------------------------------------------------------------------------------------
01999         !interior
02000         !bc=1
02001         myBC=1
02002 
02003         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
02004         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
02005 
02006         do myJ=1,maxJ
02007                 
02008                 j=jVec(myJ)
02009                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
02010                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
02011 
02012                 RHS=0
02013                 do myI=1,maxI
02014 
02015                         i=iVec(myI)
02016                         
02017                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
02018                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
02019                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
02020 
02021                         !unk1
02022                         t=1
02023                         do myT=1,endIndex(t)
02024 
02025 
02026                                 TC=x(IM(i,j,t,absolute))
02027                                 TW=x(IM(i-1,j,t,absolute))
02028                                 TE=x(IM(i+1,j,t,absolute))
02029                                 TS=x(IM(i,j-1,t,absolute))
02030                                 TN=x(IM(i,j+1,t,absolute))
02031 
02032                                 d2TdX2=(TE+TW-2.d0*TC)/(dx**2.d0)
02033 
02034                                 d2TdY2=(TN+TS-2.d0*TC)/(dy**2.d0)
02035 
02036                                 RHS(t)=d2TdX2+d2TdY2+sourceT
02037 
02038                                 !RHS(t)=TC
02039 
02040                                 !print*,TC,i,j
02041 
02042                         end do
02043 
02044                         !unk1
02045                         t=2
02046                         do myT=1,endIndex(t)
02047 
02048 
02049                                 UC=x(IM(i,j,t,absolute))
02050                                 UW=x(IM(i-1,j,t,absolute))
02051                                 UE=x(IM(i+1,j,t,absolute))
02052                                 US=x(IM(i,j-1,t,absolute))
02053                                 UN=x(IM(i,j+1,t,absolute))
02054 
02055                                 d2UdX2=(UE+UW-2.d0*UC)/(dx**2.d0)
02056 
02057                                 d2UdY2=(UN+US-2.d0*UC)/(dy**2.d0)
02058 
02059                                 RHS(t)=d2UdX2+d2UdY2+sourceU
02060 
02061                                 !RHS(t)=TC
02062 
02063                                 !print*,TC,i,j
02064 
02065                         end do
02066 
02067 
02068 
02069                         do myT=1,maxT
02070 
02071                                 t=tVec(myT)
02072 
02073                                 !print*,IM(i,j,t,relative),rhs(t)
02074 
02075                                 discrete(IM(i,j,t,relative))=rhs(t)
02076 
02077                         end do
02078 
02079                 end do
02080 
02081 
02082         end do
02083         !---------------------------------------------------------------------------------------------
02084         !stop
02085 
02086         end do ! domain 2 loop
02087 
02088 
02089 
02090 
02091 end do ! physics loop
02092 
02093 !#####################################################################
02094 
02095 
02096 deallocate(RHS)
02097 
02098 
02099 
02100 !x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x
02101 !x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x
02102 !x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x
02103 !x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x
02104 
02105 
02106 
02107 
02108 
02109 !===================================================================
02110 discreteY=>y
02111 
02112 end function discrete
02113 
02114 
02115 
02116 
02117 
02118 
02119 
02120 
02121 
02122 
02123 
02124 
02125 end module model