demona: D:/demona/2.0/manual/doc/2-31set/problems/ldc

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) PC,PW,PE,PS,PN
00053 real(8) UC,UW,UE,US,UN
00054 real(8) VC,VW,VE,VS,VN
00055 
00056 real(8) PEE,PWW,PNN,PSS
00057 
00058 real(8) UEE,UWW
00059 real(8) UEEE,UWWW
00060 real(8) UNN,USS
00061 real(8) UNNN,USSS
00062 real(8) VEE,VWW
00063 real(8) VEEE,VWWW
00064 real(8) VNN,VSS
00065 real(8) VNNN,VSSS
00066 
00067 
00068 real(8) PCPTO
00069 
00070 real(8) UCPTO,UCRTO
00071 
00072 real(8) VCPTO,VCRTO
00073 
00074 
00075 real(8) dPdX,dPdY
00076 real(8) dUdX,dUdY
00077 real(8) dVdX,dVdY
00078 real(8) d2UdX2,d2UdY2
00079 real(8) d2VdX2,d2VdY2
00080 
00081 real(8) dPdTau
00082 
00083 real(8) dUdTau,dUdT
00084 
00085 real(8) dVdTau,dVdT
00086 
00087 
00088 
00089 real(8) UdUdX,VdUdY
00090 real(8) UdVdX,VdVdY
00091 
00092 
00093 ! finite volume discretizations
00094 real(8) dUdXe,dUdXw,dUdYn,dUdYs
00095 real(8) dVdXe,dVdXw,dVdYn,dVdYs
00096 
00097 real(8) U_fd_e,U_fd_w,U_fd_n,U_fd_s
00098 real(8) V_fd_e,V_fd_w,V_fd_n,V_fd_s
00099 
00100 real(8) U_fd,V_fd
00101 
00102 real(8) U_fc_e,U_fc_w,U_fc_n,U_fc_s
00103 real(8) V_fc_e,V_fc_w,V_fc_n,V_fc_s
00104 
00105 real(8) U_fc,V_fc
00106 
00107 real(8) Ufw,Ufe,Ufn,Ufs
00108 real(8) Vfw,Vfe,Vfn,Vfs
00109 
00110 real(8) Px,Py
00111 
00112 real(8) upw1,upw2,down,value
00113 
00114 
00115 ! MLN parameters
00116 
00117 real(8) au,Ap
00118 
00119 
00120 
00121 !real(8) dtau
00122 
00123 real(8) upu,upv
00124 
00125 real(8) myU,myV
00126 
00127 real(8) upwx
00128 real(8) upwy
00129 
00130 real(8) betaAC
00131 
00132 real(8) sourceT
00133 real(8) sourceU
00134 
00135 real(8) dx,dy
00136 
00137 real(8),pointer :: myLx,myLy
00138 
00139 
00140 real(8),pointer :: phyOrder
00141 
00142 
00143 real(8),pointer :: Re
00144 
00145 
00146 !< if x+eps*v is to be evaluated, then tempy is updated before calling discrete(1)
00147 !< such that tempy(setVec)=tempy(setVec)+lambda*deltax
00148 !< if discrete(0) is called, original y is evaluated
00149 do i=1,a
00150         discreteY=>tempy
00151 end do
00152 
00153 
00154 
00155 
00156 
00157 betaAC=10.d0
00158 
00159 !betaAC=1.d0/dTau
00160 
00161 discrete=0.d0
00162 
00163 
00164 x=>discreteY
00165 
00166 xPTO=>yPTO
00167 
00168 
00169 discrete=x(setVEc)
00170 
00171 
00172 
00173 ! In this problem, discretization for 1 and 2 is the same although they are located in different
00174 ! part of the domains of different physics
00175 ! they might be coded in a loop so double explicit definition for P1 and P2 is avoided
00176 ! this is not the case for BC 3 since it is defined explicitly
00177 ! however, if and interpolation interface is used, than the bourden is passed to the interface routine
00178 ! so that the discretizations looks the same. After then, P1 and P2 can be defined in a loop.
00179 
00180 
00181 
00182 
00183 
00184 
00185 
00186 
00187 !!!!!!!!!!!!!!
00188 !< assume that  physics-1 starts here
00189 !< also that the domain-1 of physics-1 starts here, as well
00190 !#############################################################################################
00191 !#############################################################################################
00192 !#############################################################################################
00193 !
00194 !############                                   ###
00195 !#############                                  ###
00196 !###            ###                                     ###
00197 !###            ###                                     ###
00198 !###            ###                                     ###
00199 !###            ###                                     ###
00200 !#############          ####            ###
00201 !############           ####            ###
00202 !###                                                    ###
00203 !###                                                    ###
00204 !###                                                    ###
00205 !###                                                    ###
00206 !###                                                    ###
00207 !###                                                    ###     
00208 !
00209 !#############################################################################################
00210 !#############################################################################################
00211 !#############################################################################################
00212 
00213 
00214 !modelPhysics=currentPhysics
00215 !modelDomain=currentDomain
00216 
00217 modelPhysics=1
00218 
00219 
00220 !< the idea is  mymodel%discrete%set(1)
00221 !<                              mymodel%discrete%set(2)
00222 !<                              mymodel%discrete%set(3)
00223 ! with set vec pointing to the set of the physics(..)%domain%(..)%set
00224 
00225 
00226 
00227 phyOrder=>myModel%physics(modelPhysics)%evalOrder
00228 
00229 RE=>myModel%physics(modelPhysics)%physicalParameters(1)
00230 
00231 
00232 
00233 !Re=100.d0
00234 !Re=1000.d0
00235 !Re=2000.d0
00236 !Re=3000.d0
00237 !Re=3500.d0
00238 !Re=5000.d0
00239 !Re=10000.d0
00240 !stop
00241 
00242 
00243 
00244 !print*,'dx,dy=',dx,dy
00245 !stop
00246 
00247 
00248 dof=>myModel%physics(modelPhysics)%dof
00249 
00250 
00251 
00252 
00253 allocate(RHS(dof))
00254 RHS=0.d0
00255 
00256 !#####################################################################
00257 !user supplied variables
00258 !#####################################################################
00259 
00260 do myP=1,myModel%endIndex(modelPhysics)
00261 
00262         !#############################################################################################
00263         !
00264         !########                       ##
00265         !##             ##                      ##
00266         !##             ##                      ##
00267         !##             ##      ###             ##
00268         !##             ##                      ##
00269         !##             ##                      ##
00270         !########                       ##
00271         !
00272         !#############################################################################################
00273         
00274         modelDomain=1
00275 
00276         IM=>myModel%physics(modelPhysics)%domain(modelDomain)%dofMatrix
00277 
00278         myLx=>myModel%physics(modelPhysics)%domain(modelDomain)%Lx
00279         myLy=>myModel%physics(modelPhysics)%domain(modelDomain)%Ly
00280 
00281         !dx=myLx/(myModel%physics(modelPhysics)%domain(modelDomain)%maxI-myModel%physics(modelPhysics)%domain(modelDomain)%minI)
00282         !dy=myLy/(myModel%physics(modelPhysics)%domain(modelDomain)%maxJ-myModel%physics(modelPhysics)%domain(modelDomain)%minJ)
00283         dx=myLx/(myModel%physics(modelPhysics)%domain(modelDomain)%nx-1.d0)
00284         dy=myLy/(myModel%physics(modelPhysics)%domain(modelDomain)%ny-1.d0)
00285         ! normally, especially in finite volume coordinates should be read so that face lengths and center-to-center distances
00286         ! are calculated. for this uniform finite difference grid dx,dy is enough to proceed
00287 
00288         ! caution, how is dx dy in finite volume!!! was lx for the whole domain with ghostcells 
00289         ! or just the real interior domain w/o ghost cells!!! - if so number of ghostcells connection!!
00290 
00291         !print*,dx,dy
00292 
00293 
00294 
00295         do myD=1,myModel%physics(modelPhysics)%endIndex(modelDomain)
00296 
00297 
00298         !---------------------------------------------------------------------------------------------
00299         !interior
00300         !bc=1
00301         myBC=1
00302 
00303         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
00304         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
00305 
00306         do myJ=1,maxJ
00307                 
00308                 j=jVec(myJ)
00309                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
00310                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
00311 
00312                 RHS=0
00313                 do myI=1,maxI
00314 
00315                         i=iVec(myI)
00316                         
00317                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
00318                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
00319                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
00320 
00321                         !unk1
00322                         t=1
00323                         do myT=1,endIndex(t)
00324 
00325                                 PC=x(IM(i,j,t,absolute))
00326                                 !PW=x(IM(i-1,j,t,absolute))
00327                                 !PE=x(IM(i+1,j,t,absolute))
00328                                 !PN=x(IM(i,j+1,t,absolute))
00329                                 !PS=x(IM(i,j-1,t,absolute))
00330 
00331 
00332                                 UC=x(IM(i,j,t+1,absolute))
00333                                 UW=x(IM(i-1,j,t+1,absolute))
00334                                 UE=x(IM(i+1,j,t+1,absolute))
00335                                 
00336                                 VC=x(IM(i,j,t+2,absolute))
00337                                 VS=x(IM(i,j-1,t+2,absolute))
00338                                 VN=x(IM(i,j+1,t+2,absolute))
00339 
00340                                 PCPTO=xPTO(IM(i,j,t,absolute))
00341 
00342 
00343                                 !---------
00344 
00345                                 dPdTau=(PC-PCPTO)*dx*dy/dTau
00346 
00347                                 !dPdX=(PE-PW)/(2.d0*dx)
00348                                 
00349                                 !dPdY=(PN-PS)/(2.d0*dy)
00350 
00351 
00352                                 Ufw=(Uc+Uw)/2.d0
00353                                 Ufe=(Uc+Ue)/2.d0
00354 
00355                                 Vfn=(Vc+Vn)/2.d0
00356                                 Vfs=(Vc+Vs)/2.d0
00357 
00358                                 !betaAC=(max(dsqrt(10.d0),5.d0*(Uc*Uc+Vc*Vc)))**2.d0
00359 
00360                                 RHS(t)=dPdTau/betaAC+(Ufe*dy-Ufw*dy+Vfn*dx-Vfs*dx)
00361 
00362                         end do
00363 
00364                         !unk2
00365                         t=2
00366                         do myT=1,endIndex(t)
00367 
00368                                 PW=x(IM(i-1,j,t-1,absolute))
00369                                 PE=x(IM(i+1,j,t-1,absolute))
00370 
00371                                 UC=x(IM(i,j,t,absolute))
00372                                 UW=x(IM(i-1,j,t,absolute))
00373                                 UE=x(IM(i+1,j,t,absolute))
00374                                 US=x(IM(i,j-1,t,absolute))
00375                                 UN=x(IM(i,j+1,t,absolute))
00376                                 UEE=x(IM(i+2,j,t,absolute))
00377                                 UWW=x(IM(i-2,j,t,absolute))
00378                                 UNN=x(IM(i,j+2,t,absolute))
00379                                 USS=x(IM(i,j-2,t,absolute))
00380                                 
00381                                 VC=x(IM(i,j,t+1,absolute))
00382                                 VN=x(IM(i,j+1,t+1,absolute))
00383                                 VS=x(IM(i,j-1,t+1,absolute))
00384 
00385                                 UCPTO=xPTO(IM(i,j,t,absolute))
00386 
00387                                 PCPTO=xPTO(IM(i,j,t-1,absolute))
00388 
00389 
00390                                 !---------
00391 
00392                                 dPdTau=(PC-PCPTO)*dx*dy/dTau
00393                                 
00394                                 !pseudo time
00395                                 dUdTau=(UC-UCPTO)*dx*dy/dTau
00396                                 
00397                                 !pressure gradient
00398                                 Px=(PE-PW)*dy/2.d0
00399 
00400                                 !diffusion
00401                                 dUdXw=(Uc-Uw)/dx
00402                                 dUdXe=(Ue-Uc)/dx
00403                                 dUdYs=(Uc-Us)/dy
00404                                 dUdYn=(Un-Uc)/dy
00405 
00406                                 U_fd_e=dy*dUdXe
00407                                 U_fd_w=dy*dUdXw
00408                                 U_fd_n=dx*dUdYn
00409                                 U_fd_s=dx*dUdYs
00410 
00411                                 U_fd=U_fd_e-U_fd_w+U_fd_n-U_fd_s
00412 
00413 
00414                                 !convection
00415 
00416                                 Ufw=(Uc+Uw)/2.d0
00417                                 Ufe=(Uc+Ue)/2.d0
00418                                 Ufn=(Uc+Un)/2.d0
00419                                 Ufs=(Uc+Us)/2.d0
00420 
00421                                 Vfn=(Vc+Vn)/2.d0
00422                                 Vfs=(Vc+Vs)/2.d0
00423 
00424 
00425                                 !-----------------------------------------------
00426                                 !###############################################
00427                                 !-----------------------------------------------
00428                                 !U_fc_e
00429 
00430                                 if (Ufe.gt.0) then 
00431 
00432                                         down=Ue
00433                                         upw1=Uc
00434                                         upw2=Uw
00435 
00436                                 else
00437 
00438                                         down=Uc
00439                                         upw1=Ue
00440                                         upw2=UEE
00441                                 
00442                                 end if
00443                                                                 
00444                                 value=(3.d0*down+6.d0*upw1-1.d0*upw2)/8.d0
00445 
00446                                 U_fc_e=value*Ufe
00447                                 
00448                                 U_fc_e=phyOrder*U_fc_e+(1.d0-phyOrder)*Ufe*Ufe
00449                                 !-----------------------------------------------
00450                                 !###############################################
00451                                 !-----------------------------------------------
00452                                 !U_fc_w
00453 
00454 
00455                                 if (Ufw.gt.0) then 
00456 
00457                                         down=Uc
00458                                         upw1=Uw
00459                                         upw2=UWW
00460 
00461                                 else
00462 
00463                                         down=Uw
00464                                         upw1=Uc
00465                                         upw2=Ue
00466 
00467                                 end if
00468 
00469                                 value=(3.d0*down+6.d0*upw1-1.d0*upw2)/8.d0
00470                                 
00471                                 U_fc_w=value*Ufw
00472 
00473                                 U_fc_w=phyOrder*U_fc_w+(1.d0-phyOrder)*Ufw*Ufw
00474                                 !-----------------------------------------------
00475                                 !###############################################
00476                                 !-----------------------------------------------
00477                                 !U_fc_n
00478                                 
00479 
00480                                 if (Vfn.gt.0) then 
00481 
00482                                         down=Un
00483                                         upw1=Uc
00484                                         upw2=Us
00485 
00486                                 else
00487 
00488                                         down=Uc
00489                                         upw1=Un
00490                                         upw2=UNN
00491 
00492                                 end if
00493                                 
00494 
00495                                 value=(3.d0*down+6.d0*upw1-1.d0*upw2)/8.d0
00496 
00497                                 U_fc_n=value*Vfn
00498 
00499                                 U_fc_n=phyOrder*U_fc_n+(1.d0-phyOrder)*Vfn*Ufn
00500                                 !-----------------------------------------------
00501                                 !###############################################
00502                                 !-----------------------------------------------
00503                                 !U_fc_s
00504 
00505 
00506                                 if (Vfs.gt.0) then 
00507 
00508                                         down=Uc
00509                                         upw1=Us
00510                                         upw2=USS
00511 
00512                                 else
00513 
00514                                         down=Us
00515                                         upw1=Uc
00516                                         upw2=Un
00517 
00518                                 end if
00519 
00520 
00521                                 value=(3.d0*down+6.d0*upw1-1.d0*upw2)/8.d0
00522 
00523                                 U_fc_s=value*Vfs
00524 
00525                                 U_fc_s=phyOrder*U_fc_s+(1.d0-phyOrder)*Vfs*Ufs
00526 
00527                                 !-----------------------------------------------
00528                                 !###############################################
00529                                 !-----------------------------------------------
00530 
00531                         
00532                                 U_fc=U_fc_e*dy-U_fc_w*dy+U_fc_n*dx-U_fc_s*dx
00533 
00534 
00535                                 au=1.d0
00536                                 
00537                                 !betaAC=(max(dsqrt(10.d0),5.d0*(Uc*Uc+Vc*Vc)))**2.d0
00538 
00539                                 RHS(t)=(au*UC*dPdTau/betaAC+dUdTau+U_fc+Px)*Re-U_fd
00540 
00541 
00542                         end do
00543 
00544                         !unk1
00545                         t=3
00546                         do myT=1,endIndex(t)
00547 
00548                                 PC=x(IM(i,j,t-2,absolute))
00549                                 PN=x(IM(i,j+1,t-2,absolute))
00550                                 PS=x(IM(i,j-1,t-2,absolute))
00551 
00552                                 UC=x(IM(i,j,t-1,absolute))
00553                                 UW=x(IM(i-1,j,t-1,absolute))
00554                                 UE=x(IM(i+1,j,t-1,absolute))
00555 
00556                                 VC=x(IM(i,j,t,absolute))
00557                                 VW=x(IM(i-1,j,t,absolute))
00558                                 VE=x(IM(i+1,j,t,absolute))
00559                                 VS=x(IM(i,j-1,t,absolute))
00560                                 VN=x(IM(i,j+1,t,absolute))
00561                                 VWW=x(IM(i-2,j,t,absolute))
00562                                 VEE=x(IM(i+2,j,t,absolute))
00563                                 VSS=x(IM(i,j-2,t,absolute))
00564                                 VNN=x(IM(i,j+2,t,absolute))
00565                                 
00566                                 VCPTO=xPTO(IM(i,j,t,absolute))
00567 
00568                                 PCPTO=xPTO(IM(i,j,t-2,absolute))
00569 
00570 
00571                                 !---------
00572 
00573                                 dPdTau=(PC-PCPTO)*dx*dy/dTau
00574                                 !pseudo time
00575                                 dVdTau=(VC-VCPTO)*dx*dy/dTau
00576 
00577                                 !pressure gradient
00578                                 Py=(PN-PS)*dx/2.d0
00579 
00580                                 !diffusion
00581                                 dVdXw=(Vc-Vw)/dx
00582                                 dVdXe=(Ve-Vc)/dx
00583                                 dVdYs=(Vc-Vs)/dy
00584                                 dVdYn=(Vn-Vc)/dy
00585 
00586                                 V_fd_e=dy*dVdXe
00587                                 V_fd_w=dy*dVdXw
00588                                 V_fd_n=dx*dVdYn
00589                                 V_fd_s=dx*dVdYs
00590 
00591                                 V_fd=V_fd_e-V_fd_w+V_fd_n-V_fd_s
00592 
00593                                 !convection
00594 
00595                                 Ufw=(Uc+Uw)/2.d0
00596                                 Ufe=(Uc+Ue)/2.d0
00597 
00598                                 Vfw=(Vc+Vw)/2.d0
00599                                 Vfe=(Vc+Ve)/2.d0
00600                                 Vfn=(Vc+Vn)/2.d0
00601                                 Vfs=(Vc+Vs)/2.d0
00602 
00603                                 !-----------------------------------------------
00604                                 !###############################################
00605                                 !-----------------------------------------------
00606                                 !V_fc_e
00607 
00608                                 if (Ufe.gt.0) then 
00609 
00610                                         down=Ve
00611                                         upw1=Vc
00612                                         upw2=Vw
00613 
00614 
00615                                 else
00616 
00617                                         down=Vc
00618                                         upw1=Ve
00619                                         upw2=VEE
00620 
00621                                 end if
00622 
00623 
00624                                 value=(3.d0*down+6.d0*upw1-1.d0*upw2)/8.d0
00625 
00626                                 V_fc_e=value*Ufe
00627 
00628                                 V_fc_e=phyOrder*V_fc_e+(1.d0-phyOrder)*Ufe*Vfe
00629 
00630 
00631                                 !-----------------------------------------------
00632                                 !###############################################
00633                                 !-----------------------------------------------
00634                                 !V_fc_w
00635 
00636                                 if (Ufw.gt.0) then 
00637 
00638                                         down=Vc
00639                                         upw1=Vw
00640                                         upw2=VWW
00641 
00642                                 else
00643 
00644                                         down=Vw
00645                                         upw1=Vc
00646                                         upw2=Ve
00647 
00648                                 end if
00649 
00650                                 value=(3.d0*down+6.d0*upw1-1.d0*upw2)/8.d0
00651 
00652                                 V_fc_w=value*Ufw
00653 
00654                                 V_fc_w=phyOrder*V_fc_w+(1.d0-phyOrder)*Ufw*Vfw
00655 
00656                                 !-----------------------------------------------
00657                                 !###############################################
00658                                 !-----------------------------------------------
00659                                 !V_fc_n
00660 
00661                                 if (Vfn.gt.0) then 
00662 
00663                                         down=Vn
00664                                         upw1=Vc
00665                                         upw2=Vs
00666 
00667                                 else
00668 
00669                                         down=Vc
00670                                         upw1=Vn
00671                                         upw2=VNN
00672 
00673                                 end if
00674 
00675 
00676                                 value=(3.d0*down+6.d0*upw1-1.d0*upw2)/8.d0
00677 
00678                                 V_fc_n=value*Vfn
00679 
00680                                 V_fc_n=phyOrder*V_fc_n+(1.d0-phyOrder)*Vfn*Vfn
00681 
00682                                 !-----------------------------------------------
00683                                 !###############################################
00684                                 !-----------------------------------------------
00685                                 !V_fc_s
00686 
00687                                 if (Vfs.gt.0) then 
00688 
00689                                         down=Vc
00690                                         upw1=Vs
00691                                         upw2=VSS
00692                                 else
00693 
00694                                         down=Vs
00695                                         upw1=Vc
00696                                         upw2=Vn
00697 
00698 
00699                                 end if
00700 
00701 
00702                                 value=(3.d0*down+6.d0*upw1-1.d0*upw2)/8.d0
00703 
00704                                 V_fc_s=value*Vfs
00705 
00706                                 V_fc_s=phyOrder*V_fc_s+(1.d0-phyOrder)*Vfs*Vfs
00707 
00708 
00709 
00710                                 V_fc=V_fc_e*dy-V_fc_w*dy+V_fc_n*dx-V_fc_s*dx
00711 
00712 
00713                                 au=1.d0
00714 
00715                                 !betaAC=(max(dsqrt(10.d0),5.d0*(Uc*Uc+Vc*Vc)))**2.d0
00716 
00717                                 RHS(t)=(au*VC*dPdTau/betaAC+dVdTau+V_fc+Py)*Re-V_fd
00718 
00719 
00720                         end do
00721 
00722 
00723 
00724                         do myT=1,maxT
00725 
00726                                 t=tVec(myT)
00727 
00728                                 !print*,IM(i,j,t,relative),rhs(t)
00729 
00730                                 discrete(IM(i,j,t,relative))=rhs(t)
00731 
00732                         end do
00733 
00734                 end do
00735 
00736 
00737         end do
00738         !---------------------------------------------------------------------------------------------
00739         !boundary
00740         !---------------------------------------------------------------------------------------------
00741         !bottom
00742         !bc=2
00743         myBC=2 
00744 
00745         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
00746         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
00747 
00748 
00749         do myJ=1,maxJ
00750                 
00751                 j=jVec(myJ)
00752                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
00753                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
00754 
00755                 RHS=0
00756                 do myI=1,maxI
00757 
00758                         i=iVec(myI)
00759                         
00760                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
00761                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
00762                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
00763 
00764                         !unk1
00765                         t=1
00766                         do myT=1,endIndex(t)
00767 
00768                                 ! center values
00769                                 PC=x(IM(i,j,t,absolute))
00770                                 PN=x(IM(i,j+1,t,absolute))
00771 
00772                                 RHS(t)=(PN-PC)/dy
00773 
00774                         end do
00775 
00776                         !unk2
00777                         t=2
00778                         do myT=1,endIndex(t)
00779 
00780                                 ! center values
00781                                 UC=x(IM(i,j,t,absolute))
00782                                 UN=x(IM(i,j+1,t,absolute))
00783 
00784                                 !print*,TC,i,j
00785 
00786                                 RHS(t)=(UN+UC)/2.d0-0.d0
00787 
00788                         end do
00789 
00790                         !unk3
00791                         t=3
00792                         do myT=1,endIndex(t)
00793 
00794                                 ! center values
00795                                 VC=x(IM(i,j,t,absolute))
00796                                 VN=x(IM(i,j+1,t,absolute))
00797 
00798                                 !print*,TC,i,j
00799 
00800                                 RHS(t)=(VN+VC)/2.d0-0.d0
00801 
00802                         end do
00803 
00804 
00805                         do myT=1,maxT
00806                                 
00807                                 t=tVec(myT)
00808 
00809                                 !print*,IM(i,j,t,relative)
00810 
00811                                 discrete(IM(i,j,t,relative))=rhs(t)
00812 
00813 
00814                         end do
00815 
00816                 end do
00817 
00818 
00819         end do
00820         !---------------------------------------------------------------------------------------------
00821         !stop
00822         !---------------------------------------------------------------------------------------------
00823         !top
00824         !bc=3
00825         myBC=3 
00826 
00827         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
00828         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
00829 
00830 
00831         do myJ=1,maxJ
00832                 
00833                 j=jVec(myJ)
00834                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
00835                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
00836 
00837                 RHS=0
00838                 do myI=1,maxI
00839 
00840                         i=iVec(myI)
00841                         
00842                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
00843                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
00844                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
00845 
00846                         !unk1
00847                         t=1
00848                         do myT=1,endIndex(t)
00849 
00850                                 ! center values
00851                                 PC=x(IM(i,j,t,absolute))
00852                                 PS=x(IM(i,j-1,t,absolute))
00853 
00854                                 RHS(t)=(PC-PS)/dy
00855 
00856                         end do
00857 
00858 
00859                         !unk2
00860                         t=2
00861                         do myT=1,endIndex(t)
00862 
00863                                 ! center values
00864                                 UC=x(IM(i,j,t,absolute))
00865                                 US=x(IM(i,j-1,t,absolute))
00866 
00867                                 !print*,TC,i,j
00868 
00869                                 RHS(t)=(UC+US)/2.d0-1.d0
00870 
00871                         end do
00872 
00873                         !unk3
00874                         t=3
00875                         do myT=1,endIndex(t)
00876 
00877                                 ! center values
00878                                 VC=x(IM(i,j,t,absolute))
00879                                 VS=x(IM(i,j-1,t,absolute))
00880 
00881                                 !print*,TC,i,j
00882 
00883                                 RHS(t)=(VC+VS)/2.d0-0.d0
00884 
00885                         end do
00886 
00887 
00888                         do myT=1,maxT
00889                                 
00890                                 t=tVec(myT)
00891 
00892                                 !print*,IM(i,j,t,relative)
00893 
00894                                 discrete(IM(i,j,t,relative))=rhs(t)
00895 
00896 
00897                         end do
00898 
00899                 end do
00900 
00901 
00902         end do
00903         !---------------------------------------------------------------------------------------------
00904         !stop
00905         !---------------------------------------------------------------------------------------------
00906         !left
00907         !bc=4
00908         myBC=4   
00909 
00910         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
00911         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
00912 
00913 
00914         do myJ=1,maxJ
00915                 
00916                 j=jVec(myJ)
00917                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
00918                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
00919 
00920                 RHS=0
00921                 do myI=1,maxI
00922 
00923                         i=iVec(myI)
00924                         
00925                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
00926                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
00927                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
00928 
00929 
00930 
00931                         !unk1
00932                         t=1
00933                         do myT=1,endIndex(t)
00934 
00935                                 ! center values
00936                                 PC=x(IM(i,j,t,absolute))
00937                                 PE=x(IM(i+1,j,t,absolute))
00938                                 
00939                                 RHS(t)=(PE-PC)/dx
00940 
00941                         end do
00942 
00943                         !unk2
00944                         t=2
00945                         do myT=1,endIndex(t)
00946 
00947                                 ! center values
00948                                 UC=x(IM(i,j,t,absolute)) !< or TC=x(IM(i,j,*1*,absolute))
00949                                 UE=x(IM(i+1,j,t,absolute)) !< or TC=x(IM(i,j,*1*,absolute))
00950 
00951                                 RHS(t)=(UE+UC)/2.d0-0.d0
00952 
00953                         end do
00954 
00955                         !unk3
00956                         t=3
00957                         do myT=1,endIndex(t)
00958 
00959                                 ! center values
00960                                 VC=x(IM(i,j,t,absolute))
00961                                 VE=x(IM(i+1,j,t,absolute))
00962 
00963                                 RHS(t)=(VE+VC)/2.d0-0.d0
00964 
00965                         end do
00966 
00967 
00968                         do myT=1,maxT
00969                                 
00970                                 t=tVec(myT)
00971 
00972                                 !print*,IM(i,j,t,relative)
00973 
00974                                 discrete(IM(i,j,t,relative))=rhs(t)
00975 
00976                         end do
00977 
00978                 end do
00979 
00980 
00981 
00982         end do
00983         !---------------------------------------------------------------------------------------------
00984         !stop
00985         !---------------------------------------------------------------------------------------------
00986         !right
00987         !bc=5
00988         myBC=5 
00989 
00990         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
00991         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
00992 
00993 
00994         do myJ=1,maxJ
00995                 
00996                 j=jVec(myJ)
00997                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
00998                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
00999 
01000                 RHS=0
01001                 do myI=1,maxI
01002 
01003                         i=iVec(myI)
01004                         
01005                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
01006                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
01007                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
01008 
01009                         !unk1
01010                         t=1
01011                         do myT=1,endIndex(t)
01012 
01013                                 ! center values
01014                                 PC=x(IM(i,j,t,absolute))
01015                                 PW=x(IM(i-1,j,t,absolute))
01016 
01017                                 !print*,TC,i,j
01018 
01019                                 RHS(t)=(PC-PW)/dx
01020 
01021                         end do
01022 
01023                         !unk2
01024                         t=2
01025                         do myT=1,endIndex(t)
01026 
01027                                 ! center values
01028                                 UC=x(IM(i,j,t,absolute))
01029                                 UW=x(IM(i-1,j,t,absolute))
01030 
01031                                 !print*,TC,i,j
01032 
01033                                 RHS(t)=(UC+UW)/2.d0-0.d0
01034 
01035                         end do
01036 
01037                         !unk3
01038                         t=3
01039                         do myT=1,endIndex(t)
01040 
01041                                 ! center values
01042                                 VC=x(IM(i,j,t,absolute))
01043                                 VW=x(IM(i-1,j,t,absolute))
01044 
01045                                 !print*,TC,i,j
01046 
01047                                 RHS(t)=(VC+VW)/2.d0-0.d0
01048 
01049                         end do
01050 
01051 
01052                         do myT=1,maxT
01053                                 
01054                                 t=tVec(myT)
01055 
01056                                 !print*,IM(i,j,t,relative)
01057 
01058                                 discrete(IM(i,j,t,relative))=rhs(t)
01059 
01060 
01061                         end do
01062 
01063                 end do
01064 
01065 
01066         end do
01067         !---------------------------------------------------------------------------------------------
01068         !stop
01069         !---------------------------------------------------------------------------------------------
01070         !bottom - ghost cell 2
01071         !bc=6
01072         myBC=6 
01073 
01074         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
01075         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
01076 
01077 
01078         do myJ=1,maxJ
01079                 
01080                 j=jVec(myJ)
01081                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
01082                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
01083 
01084                 RHS=0
01085                 do myI=1,maxI
01086 
01087                         i=iVec(myI)
01088                         
01089                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
01090                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
01091                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
01092 
01093                         !unk1
01094                         t=1
01095                         do myT=1,endIndex(t)
01096 
01097                                 ! center values
01098                                 PC=x(IM(i,j,t,absolute))
01099                                 PN=x(IM(i,j+1,t,absolute))
01100 
01101                                 RHS(t)=(PN-PC)/dy
01102 
01103                         end do
01104 
01105                         !unk2
01106                         t=2
01107                         do myT=1,endIndex(t)
01108 
01109                                 ! center values
01110                                 UC=x(IM(i,j,t,absolute))
01111                                 UN=x(IM(i,j+1,t,absolute))
01112                                 UNNN=x(IM(i,j+3,t,absolute))
01113 
01114                                 RHS(t)=phyOrder*((UNNN+UC)/2.d0-0.d0) + (1.d0-phyOrder)*(UC-3.d0*UN+2.d0*0.d0)
01115 
01116                         end do
01117 
01118                         !unk3
01119                         t=3
01120                         do myT=1,endIndex(t)
01121 
01122                                 ! center values
01123                                 VC=x(IM(i,j,t,absolute))
01124                                 VN=x(IM(i,j+1,t,absolute))
01125                                 VNNN=x(IM(i,j+3,t,absolute))
01126 
01127                                 RHS(t)=phyOrder*((VNNN+VC)/2.d0-0.d0) + (1.d0-phyOrder)*(VC-3.d0*VN+2.d0*0.d0)
01128 
01129                         end do
01130 
01131 
01132                         do myT=1,maxT
01133                                 
01134                                 t=tVec(myT)
01135 
01136                                 !print*,IM(i,j,t,relative)
01137 
01138                                 discrete(IM(i,j,t,relative))=rhs(t)
01139 
01140 
01141                         end do
01142 
01143                 end do
01144 
01145 
01146         end do
01147         !---------------------------------------------------------------------------------------------
01148         !stop
01149         !---------------------------------------------------------------------------------------------
01150         !top - ghost cell 2
01151         !bc=7
01152         myBC=7 
01153 
01154         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
01155         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
01156 
01157 
01158         do myJ=1,maxJ
01159                 
01160                 j=jVec(myJ)
01161                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
01162                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
01163 
01164                 RHS=0
01165                 do myI=1,maxI
01166 
01167                         i=iVec(myI)
01168                         
01169                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
01170                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
01171                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
01172 
01173                         !unk1
01174                         t=1
01175                         do myT=1,endIndex(t)
01176 
01177                                 ! center values
01178                                 PC=x(IM(i,j,t,absolute))
01179                                 PS=x(IM(i,j-1,t,absolute))
01180 
01181                                 RHS(t)=(PC-PS)/dy
01182 
01183                         end do
01184 
01185 
01186                         !unk2
01187                         t=2
01188                         do myT=1,endIndex(t)
01189 
01190                                 ! center values
01191                                 UC=x(IM(i,j,t,absolute))
01192                                 US=x(IM(i,j-1,t,absolute))
01193                                 USSS=x(IM(i,j-3,t,absolute))
01194 
01195                                 RHS(t)=phyOrder*((USSS+UC)/2.d0-1.d0) + (1.d0-phyOrder)*(UC-3.d0*US+2.d0*1.d0)
01196 
01197                         end do
01198 
01199                         !unk3
01200                         t=3
01201                         do myT=1,endIndex(t)
01202 
01203                                 ! center values
01204                                 VC=x(IM(i,j,t,absolute))
01205                                 VS=x(IM(i,j-1,t,absolute))
01206                                 VSSS=x(IM(i,j-3,t,absolute))
01207 
01208                                 !print*,TC,i,j
01209 
01210                                 RHS(t)=phyOrder*((VSSS+VC)/2.d0-0.d0) + (1.d0-phyOrder)*(VC-3.d0*VS+2.d0*0.d0)
01211 
01212                         end do
01213 
01214 
01215                         do myT=1,maxT
01216                                 
01217                                 t=tVec(myT)
01218 
01219                                 !print*,IM(i,j,t,relative)
01220 
01221                                 discrete(IM(i,j,t,relative))=rhs(t)
01222 
01223 
01224                         end do
01225 
01226                 end do
01227 
01228 
01229         end do
01230         !---------------------------------------------------------------------------------------------
01231         !stop
01232         !---------------------------------------------------------------------------------------------
01233         !left - ghost cell 2
01234         !bc=8
01235         myBC=8   
01236 
01237         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
01238         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
01239 
01240 
01241         do myJ=1,maxJ
01242                 
01243                 j=jVec(myJ)
01244                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
01245                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
01246 
01247                 RHS=0
01248                 do myI=1,maxI
01249 
01250                         i=iVec(myI)
01251                         
01252                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
01253                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
01254                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
01255 
01256 
01257 
01258                         !unk1
01259                         t=1
01260                         do myT=1,endIndex(t)
01261 
01262                                 ! center values
01263                                 PC=x(IM(i,j,t,absolute))
01264                                 PE=x(IM(i+1,j,t,absolute))
01265                                 
01266                                 RHS(t)=(PE-PC)/dx
01267 
01268                         end do
01269 
01270                         !unk2
01271                         t=2
01272                         do myT=1,endIndex(t)
01273 
01274                                 ! center values
01275                                 UC=x(IM(i,j,t,absolute)) !< or TC=x(IM(i,j,*1*,absolute))
01276                                 UE=x(IM(i+1,j,t,absolute)) !< or TC=x(IM(i,j,*1*,absolute))
01277                                 UEEE=x(IM(i+3,j,t,absolute)) !< or TC=x(IM(i,j,*1*,absolute))
01278 
01279                                 RHS(t)=phyOrder*((UEEE+UC)/2.d0-0.d0) + (1.d0-phyOrder)*(UC-3.d0*UE+2.d0*0.d0)
01280 
01281                         end do
01282 
01283                         !unk3
01284                         t=3
01285                         do myT=1,endIndex(t)
01286 
01287                                 ! center values
01288                                 VC=x(IM(i,j,t,absolute))
01289                                 VE=x(IM(i+1,j,t,absolute))
01290                                 VEEE=x(IM(i+3,j,t,absolute))
01291 
01292                                 RHS(t)=phyOrder*((VEEE+VC)/2.d0-0.d0) + (1.d0-phyOrder)*(VC-3.d0*VE+2.d0*0.d0)
01293 
01294                         end do
01295 
01296 
01297                         do myT=1,maxT
01298                                 
01299                                 t=tVec(myT)
01300 
01301                                 !print*,IM(i,j,t,relative)
01302 
01303                                 discrete(IM(i,j,t,relative))=rhs(t)
01304 
01305                         end do
01306 
01307                 end do
01308 
01309 
01310 
01311         end do
01312         !---------------------------------------------------------------------------------------------
01313         !stop
01314         !---------------------------------------------------------------------------------------------
01315         !right - ghost cell 2
01316         !bc=9
01317         myBC=9 
01318 
01319         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
01320         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
01321 
01322 
01323         do myJ=1,maxJ
01324                 
01325                 j=jVec(myJ)
01326                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
01327                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
01328 
01329                 RHS=0
01330                 do myI=1,maxI
01331 
01332                         i=iVec(myI)
01333                         
01334                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
01335                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
01336                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
01337 
01338                         !unk1
01339                         t=1
01340                         do myT=1,endIndex(t)
01341 
01342                                 ! center values
01343                                 PC=x(IM(i,j,t,absolute))
01344                                 PW=x(IM(i-1,j,t,absolute))
01345 
01346                                 !print*,TC,i,j
01347 
01348                                 RHS(t)=(PC-PW)/dx
01349 
01350                         end do
01351 
01352                         !unk2
01353                         t=2
01354                         do myT=1,endIndex(t)
01355 
01356                                 ! center values
01357                                 UC=x(IM(i,j,t,absolute))
01358                                 UW=x(IM(i-1,j,t,absolute))
01359                                 UWWW=x(IM(i-3,j,t,absolute))
01360 
01361                                 !print*,TC,i,j
01362 
01363                                 RHS(t)=phyOrder*((UWWW+UC)/2.d0-0.d0) + (1.d0-phyOrder)*(UC-3.d0*UW+2.d0*0.d0)
01364 
01365                         end do
01366 
01367                         !unk3
01368                         t=3
01369                         do myT=1,endIndex(t)
01370 
01371                                 ! center values
01372                                 VC=x(IM(i,j,t,absolute))
01373                                 VW=x(IM(i-1,j,t,absolute))
01374                                 VWWW=x(IM(i-3,j,t,absolute))
01375 
01376                                 !print*,TC,i,j
01377 
01378                                 RHS(t)=phyOrder*((VWWW+VC)/2.d0-0.d0) + (1.d0-phyOrder)*(VC-3.d0*VW+2.d0*0.d0)
01379 
01380                         end do
01381 
01382 
01383                         do myT=1,maxT
01384                                 
01385                                 t=tVec(myT)
01386 
01387                                 !print*,IM(i,j,t,relative)
01388 
01389                                 discrete(IM(i,j,t,relative))=rhs(t)
01390 
01391 
01392                         end do
01393 
01394                 end do
01395 
01396 
01397         end do
01398         !---------------------------------------------------------------------------------------------
01399 
01400         end do ! domain 1 loop
01401 
01402 
01403 
01404 
01405 end do ! physics loop
01406 
01407 !#####################################################################
01408 
01409 
01410 deallocate(RHS)
01411 
01412 
01413 
01414 
01415 
01416 
01417 !===================================================================
01418 discreteY=>y
01419 
01420 end function discrete
01421 
01422 
01423 
01424 
01425 
01426 
01427 
01428 
01429 
01430 
01431 
01432 
01433 end module model
D:/demona/2.0/manual/doc/2-31set/problems/ldc_fv/1/modelMLN1.f90
