最近接到一个任务，就是解决FVCOM中对流扩散计算不守衡问题。导师认为是其求解时候水平和垂向计算分开求解所导致的，目前我也没搞清到底有什么问题，反正就是让把SUNTANS的对流扩散计算挪到FVCOM中，下面就把 SUNTANS 和 FVCOM 数值求解的过程贴出来，备忘
SUNTANS模型求解过程
SUNTANS模型手册：http://web.stanford.edu/group/suntans/cgi-bin/documentation/user_guide/user_guide.html
介绍文献：《An unstructured-grid, finite-volume, nonhydrostatic, parallel coastal ocean simulator》
代码所谓研究讨论之用这里只公布部分：
1 /*
2 * File: scalars.c
3 * Author: Oliver B. Fringer
4 * Institution: Stanford University
5 * ----------------------------------------
6 * This file contains the scalar transport function.
7 *
8 * Copyright (C) 2005-2006 The Board of Trustees of the Leland Stanford Junior
9 * University. All Rights Reserved.
10 *
11 */
12 #include "scalars.h"
13 #include "util.h"
14 #include "tvd.h"
15 #include "initialization.h"
16
17 #define SMALL_CONSISTENCY 1e-5
18
19 REAL smin_value, smax_value;
20
21 /*
22 * Function: UpdateScalars
23 * Usage: UpdateScalars(grid,phys,prop,wnew,scalar,Cn,kappa,kappaH,kappa_tv,theta);
24 * ---------------------------------------------------------------------------
25 * Update the scalar quantity stored in the array denoted by scal using the
26 * theta method for vertical advection and vertical diffusion and Adams-Bashforth
27 * for horizontal advection and diffusion.
28 *
29 * Cn must store the AB terms from time step n-1 for this scalar
30 * kappa denotes the vertical scalar diffusivity
31 * kappaH denotes the horizontal scalar diffusivity
32 * kappa_tv denotes the vertical turbulent scalar diffusivity
33 *
34 */
35 void UpdateScalars(gridT *grid, physT *phys, propT *prop, REAL **wnew, REAL **scal, REAL **boundary_scal, REAL **Cn,
36 REAL kappa, REAL kappaH, REAL **kappa_tv, REAL theta,
37 REAL **src1, REAL **src2, REAL *Ftop, REAL *Fbot, int alpha_top, int alpha_bot,
38 MPI_Comm comm, int myproc, int checkflag, int TVDscheme)
39 {
40 int i, iptr, j, jptr, ib, k, nf, ktop;
41 int Nc=grid->Nc, normal, nc1, nc2, ne;
42 REAL df, dg, Ac, dt=prop->dt, fab, *a, *b, *c, *d, *ap, *am, *bd, *uflux, dznew, mass, *sp, *temp;
43 REAL smin, smax, div_local, div_da;
44 int k1, k2, kmin, imin, kmax, imax, mincount, maxcount, allmincount, allmaxcount, flag;
45
46 prop->TVD = TVDscheme;
47 // These are used mostly debugging to turn on/off vertical and horizontal TVD.
48 prop->horiTVD = 1;
49 prop->vertTVD = 1;
50
51 ap = phys->ap;
52 am = phys->am;
53 bd = phys->bp;
54 temp = phys->bm;
55 a = phys->a;
56 b = phys->b;
57 c = phys->c;
58 d = phys->d;
59
60 // Never use AB2
61 if(1) {
62 fab=1;
63 for(i=0;i<grid->Nc;i++)
64 for(k=0;k<grid->Nk[i];k++)
65 Cn[i][k]=0;
66 } else
67 fab=1.5;
68
69 for(i=0;i<Nc;i++)
70 for(k=0;k<grid->Nk[i];k++)
71 phys->stmp[i][k]=scal[i][k];
72
73 // Add on boundary fluxes, using stmp2 as the temporary storage
74 // variable
75 //for(iptr=grid->celldist[0];iptr<grid->celldist[1];iptr++) {
76 for(iptr=grid->celldist[0];iptr<grid->celldist[2];iptr++) {
77 i = grid->cellp[iptr];
78
79 for(k=grid->ctop[i];k<grid->Nk[i];k++)
80 phys->stmp2[i][k]=0;
81 }
82
83 if(boundary_scal) {
84 for(jptr=grid->edgedist[2];jptr<grid->edgedist[5];jptr++) {
85 j = grid->edgep[jptr];
86
87 ib = grid->grad[2*j];
88
89 // Set the value of stmp2 adjacent to the boundary to the value of the boundary.
90 // This will be used to add the boundary flux when stmp2 is used again below.
91 for(k=grid->ctop[ib];k<grid->Nk[ib];k++)
92 phys->stmp2[ib][k]=boundary_scal[jptr-grid->edgedist[2]][k];
93 }
94 }
95
96 // Compute the scalar on the vertical faces (for horiz. advection)
97
98 if(prop->TVD && prop->horiTVD)
99 HorizontalFaceScalars(grid,phys,prop,scal,boundary_scal,prop->TVD,comm,myproc);
100
101 //for(iptr=grid->celldist[0];iptr<grid->celldist[1];iptr++) {
102 for(iptr=grid->celldist[0];iptr<grid->celldist[2];iptr++) {
103 i = grid->cellp[iptr];
104 Ac = grid->Ac[i];
105
106 if(grid->ctop[i]>=grid->ctopold[i]) {
107 ktop=grid->ctop[i];
108 dznew=grid->dzz[i][ktop];
109 } else {
110 ktop=grid->ctopold[i];
111 dznew=0;
112 for(k=grid->ctop[i];k<=grid->ctopold[i];k++)
113 dznew+=grid->dzz[i][k];
114 }
115
116 // These are the advective components of the tridiagonal
117 // at the new time step.
118 if(!(prop->TVD && prop->vertTVD))
119 for(k=0;k<grid->Nk[i]+1;k++) {
120 ap[k] = 0.5*(wnew[i][k]+fabs(wnew[i][k]));
121 am[k] = 0.5*(wnew[i][k]-fabs(wnew[i][k]));
122 }
123 else // Compute the ap/am for TVD schemes
124 GetApAm(ap,am,phys->wp,phys->wm,phys->Cp,phys->Cm,phys->rp,phys->rm,
125 wnew,grid->dzz,scal,i,grid->Nk[i],ktop,prop->dt,prop->TVD);
126
127 for(k=ktop+1;k<grid->Nk[i];k++) {
128 a[k-ktop]=theta*dt*am[k];
129 b[k-ktop]=grid->dzz[i][k]+theta*dt*(ap[k]-am[k+1]);
130 c[k-ktop]=-theta*dt*ap[k+1];
131 }
132
133 // Top cell advection
134 a[0]=0;
135 b[0]=dznew-theta*dt*am[ktop+1];
136 c[0]=-theta*dt*ap[ktop+1];
137
138 // Bottom cell no-flux boundary condition for advection
139 b[(grid->Nk[i]-1)-ktop]+=c[(grid->Nk[i]-1)-ktop];
140
141 // Implicit vertical diffusion terms
142 for(k=ktop+1;k<grid->Nk[i];k++)
143 bd[k]=(2.0*kappa+kappa_tv[i][k-1]+kappa_tv[i][k])/
144 (grid->dzz[i][k-1]+grid->dzz[i][k]);
145
146 for(k=ktop+1;k<grid->Nk[i]-1;k++) {
147 a[k-ktop]-=theta*dt*bd[k];
148 b[k-ktop]+=theta*dt*(bd[k]+bd[k+1]);
149 c[k-ktop]-=theta*dt*bd[k+1];
150 }
151 if(src1)
152 for(k=ktop;k<grid->Nk[i];k++)
153 b[k-ktop]+=grid->dzz[i][k]*src1[i][k]*theta*dt;
154
155 // Diffusive fluxes only when more than 1 layer
156 if(ktop<grid->Nk[i]-1) {
157 // Top cell diffusion
158 b[0]+=theta*dt*(bd[ktop+1]+2*alpha_top*bd[ktop+1]);
159 c[0]-=theta*dt*bd[ktop+1];
160
161 // Bottom cell diffusion
162 a[(grid->Nk[i]-1)-ktop]-=theta*dt*bd[grid->Nk[i]-1];
163 b[(grid->Nk[i]-1)-ktop]+=theta*dt*(bd[grid->Nk[i]-1]+2*alpha_bot*bd[grid->Nk[i]-1]);
164 }
165
166 // Explicit part into source term d[]
167 for(k=ktop+1;k<grid->Nk[i];k++)
168 d[k-ktop]=grid->dzzold[i][k]*phys->stmp[i][k];
169 if(src1)
170 for(k=ktop+1;k<grid->Nk[i];k++)
171 d[k-ktop]-=src1[i][k]*(1-theta)*dt*grid->dzzold[i][k]*phys->stmp[i][k];
172
173 d[0]=0;
174 if(grid->ctopold[i]<=grid->ctop[i]) {
175 for(k=grid->ctopold[i];k<=grid->ctop[i];k++)
176 d[0]+=grid->dzzold[i][k]*phys->stmp[i][k];
177 if(src1)
178 for(k=grid->ctopold[i];k<=grid->ctop[i];k++)
179 d[0]-=src1[i][k]*(1-theta)*dt*grid->dzzold[i][k]*phys->stmp[i][k];
180 } else {
181 d[0]=grid->dzzold[i][ktop]*phys->stmp[i][ktop];
182 if(src1)
183 d[0]-=src1[i][ktop]*(1-theta)*dt*grid->dzzold[i][ktop]*phys->stmp[i][k];
184 }
185
186 // These are the advective components of the tridiagonal
187 // that use the new velocity
188 if(!(prop->TVD && prop->vertTVD))
189 for(k=0;k<grid->Nk[i]+1;k++) {
190 ap[k] = 0.5*(phys->wtmp2[i][k]+fabs(phys->wtmp2[i][k]));
191 am[k] = 0.5*(phys->wtmp2[i][k]-fabs(phys->wtmp2[i][k]));
192 }
193 else // Compute the ap/am for TVD schemes
194 GetApAm(ap,am,phys->wp,phys->wm,phys->Cp,phys->Cm,phys->rp,phys->rm,
195 phys->wtmp2,grid->dzzold,phys->stmp,i,grid->Nk[i],ktop,prop->dt,prop->TVD);
196
197 // Explicit advection and diffusion
198 for(k=ktop+1;k<grid->Nk[i]-1;k++)
199 d[k-ktop]-=(1-theta)*dt*(am[k]*phys->stmp[i][k-1]+
200 (ap[k]-am[k+1])*phys->stmp[i][k]-
201 ap[k+1]*phys->stmp[i][k+1])-
202 (1-theta)*dt*(bd[k]*phys->stmp[i][k-1]
203 -(bd[k]+bd[k+1])*phys->stmp[i][k]
204 +bd[k+1]*phys->stmp[i][k+1]);
205
206 if(ktop<grid->Nk[i]-1) {
207 //Flux through bottom of top cell
208 k=ktop;
209 d[0]=d[0]-(1-theta)*dt*(-am[k+1]*phys->stmp[i][k]-
210 ap[k+1]*phys->stmp[i][k+1])+
211 (1-theta)*dt*(-(2*alpha_top*bd[k+1]+bd[k+1])*phys->stmp[i][k]+
212 bd[k+1]*phys->stmp[i][k+1]);
213 if(Ftop) d[0]+=dt*(1-alpha_top+2*alpha_top*bd[k+1])*Ftop[i];
214
215 // Through top of bottom cell
216 k=grid->Nk[i]-1;
217 d[k-ktop]-=(1-theta)*dt*(am[k]*phys->stmp[i][k-1]+
218 ap[k]*phys->stmp[i][k])-
219 (1-theta)*dt*(bd[k]*phys->stmp[i][k-1]-
220 (bd[k]+2*alpha_bot*bd[k])*phys->stmp[i][k]);
221 if(Fbot) d[k-ktop]+=dt*(-1+alpha_bot+2*alpha_bot*bd[k])*Fbot[i];
222 }
223
224 // First add on the source term from the previous time step.
225 if(grid->ctop[i]<=grid->ctopold[i]) {
226 for(k=grid->ctop[i];k<=grid->ctopold[i];k++)
227 d[0]+=(1-fab)*Cn[i][grid->ctopold[i]]/(1+fabs(grid->ctop[i]-grid->ctopold[i]));
228 for(k=grid->ctopold[i]+1;k<grid->Nk[i];k++)
229 d[k-grid->ctopold[i]]+=(1-fab)*Cn[i][k];
230 } else {
231 for(k=grid->ctopold[i];k<=grid->ctop[i];k++)
232 d[0]+=(1-fab)*Cn[i][k];
233 for(k=grid->ctop[i]+1;k<grid->Nk[i];k++)
234 d[k-grid->ctop[i]]+=(1-fab)*Cn[i][k];
235 }
236
237 for(k=0;k<grid->ctop[i];k++)
238 Cn[i][k]=0;
239
240 if(src2)
241 for(k=grid->ctop[i];k<grid->Nk[i];k++)
242 Cn[i][k-ktop]=dt*src2[i][k]*grid->dzzold[i][k];
243 else
244 for(k=grid->ctop[i];k<grid->Nk[i];k++)
245 Cn[i][k]=0;
246
247 // Now create the source term for the current time step
248 for(k=0;k<grid->Nk[i];k++)
249 ap[k]=0;
250
251 for(nf=0;nf<grid->nfaces[i];nf++) {
252 ne = grid->face[i*grid->maxfaces+nf];
253 normal = grid->normal[i*grid->maxfaces+nf];
254 df = grid->df[ne];
255 dg = grid->dg[ne];
256 nc1 = grid->grad[2*ne];
257 nc2 = grid->grad[2*ne+1];
258 if(nc1==-1) nc1=nc2;
259 if(nc2==-1) {
260 nc2=nc1;
261 if(boundary_scal && (grid->mark[ne]==2 || grid->mark[ne]==3))
262 sp=phys->stmp2[nc1];
263 else
264 sp=phys->stmp[nc1];
265 } else
266 sp=phys->stmp[nc2];
267
268 if(!(prop->TVD && prop->horiTVD)) {
269 for(k=0;k<grid->Nke[ne];k++)
270 temp[k]=UpWind(phys->utmp2[ne][k],
271 phys->stmp[nc1][k],
272 sp[k]);
273 } else {
274 for(k=0;k<grid->Nke[ne];k++)
275 if(phys->utmp2[ne][k]>0)
276 temp[k]=phys->SfHp[ne][k];
277 else
278 temp[k]=phys->SfHm[ne][k];
279 }
280
281 for(k=0;k<grid->Nke[ne];k++)
282 ap[k] += dt*df*normal/Ac*(theta*phys->u[ne][k]+(1-theta)*phys->utmp2[ne][k])
283 *temp[k]*grid->dzf[ne][k];
284 }
285
286 for(k=ktop+1;k<grid->Nk[i];k++)
287 Cn[i][k-ktop]-=ap[k];
288
289 for(k=0;k<=ktop;k++)
290 Cn[i][0]-=ap[k];
291
292 // Add on the source from the current time step to the rhs.
293 for(k=0;k<grid->Nk[i]-ktop;k++)
294 d[k]+=fab*Cn[i][k];
295
296 // Add on the volume correction if h was < -d
297 /*
298 if(grid->ctop[i]==grid->Nk[i]-1)
299 d[grid->Nk[i]-ktop-1]+=phys->hcorr[i]*phys->stmp[i][grid->ctop[i]];
300 */
301
302 for(k=ktop;k<grid->Nk[i];k++)
303 ap[k]=Cn[i][k-ktop];
304 for(k=0;k<=ktop;k++)
305 Cn[i][k]=0;
306 for(k=ktop+1;k<grid->Nk[i];k++)
307 Cn[i][k]=ap[k];
308 for(k=grid->ctop[i];k<=ktop;k++)
309 Cn[i][k]=ap[ktop]/(1+fabs(grid->ctop[i]-ktop));
310
311 if(grid->Nk[i]-ktop>1)
312 TriSolve(a,b,c,d,&(scal[i][ktop]),grid->Nk[i]-ktop);
313 else if(prop->n>1) {
314 if(b[0]>0 && phys->active[i])
315 scal[i][ktop]=d[0]/b[0];
316 else
317 scal[i][ktop]=0;
318 }
319
320 for(k=0;k<grid->ctop[i];k++)
321 scal[i][k]=0;
322
323 for(k=grid->ctop[i];k<grid->ctopold[i];k++)
324 scal[i][k]=scal[i][ktop];
325 }
326
327 // Code to check divergence change CHECKCONSISTENCY to 1 in suntans.h
328 if(CHECKCONSISTENCY && checkflag) {
329
330 if(prop->n==1+prop->nstart) {
331 smin=INFTY;
332 smax=-INFTY;
333 for(i=0;i<grid->Nc;i++) {
334 for(k=grid->ctop[i];k<grid->Nk[i];k++) {
335 if(phys->stmp[i][k]>smax) {
336 smax=phys->stmp[i][k];
337 imax=i;
338 kmax=k;
339 }
340 if(phys->stmp[i][k]<smin) {
341 smin=phys->stmp[i][k];
342 imin=i;
343 kmin=k;
344 }
345 }
346 }
347 MPI_Reduce(&smin,&smin_value,1,MPI_DOUBLE,MPI_MIN,0,comm);
348 MPI_Reduce(&smax,&smax_value,1,MPI_DOUBLE,MPI_MAX,0,comm);
349 MPI_Bcast(&smin_value,1,MPI_DOUBLE,0,comm);
350 MPI_Bcast(&smax_value,1,MPI_DOUBLE,0,comm);
351
352 if(myproc==0)
353 printf("Minimum scalar: %.2f, maximum: %.2f\n",smin_value,smax_value);
354 }
355
356 //for(iptr=grid->celldist[0];iptr<grid->celldist[1];iptr++) {
357 for(iptr=grid->celldist[0];iptr<grid->celldist[2];iptr++) {
358 i = grid->cellp[iptr];
359
360 flag=0;
361 for(nf=0;nf<grid->nfaces[i];nf++) {
362 if(grid->mark[grid->face[i*grid->maxfaces+nf]]==2 ||
363 grid->mark[grid->face[i*grid->maxfaces+nf]]==3) {
364 flag=1;
365 break;
366 }
367 }
368
369 if(!flag) {
370 div_da=0;
371
372 for(k=0;k<grid->Nk[i];k++) {
373 div_da+=grid->Ac[i]*(grid->dzz[i][k]-grid->dzzold[i][k])/prop->dt;
374
375 div_local=0;
376 for(nf=0;nf<grid->nfaces[i];nf++) {
377 ne=grid->face[i*grid->maxfaces+nf];
378 div_local+=(theta*phys->u[ne][k]+(1-theta)*phys->utmp2[ne][k])
379 *grid->dzf[ne][k]*grid->normal[i*grid->maxfaces+nf]*grid->df[ne];
380 }
381 div_da+=div_local;
382 div_local+=grid->Ac[i]*(theta*(wnew[i][k]-wnew[i][k+1])+
383 (1-theta)*(phys->wtmp2[i][k]-phys->wtmp2[i][k+1]));
384
385 if(k>=grid->ctop[i]) {
386 if(fabs(div_local)>SMALL_CONSISTENCY && grid->dzz[imin][0]>DRYCELLHEIGHT)
387 printf("Step: %d, proc: %d, locally-divergent at %d, %d, div=%e\n",
388 prop->n,myproc,i,k,div_local);
389 }
390 }
391 if(fabs(div_da)>SMALL_CONSISTENCY && phys->h[i]+grid->dv[i]>DRYCELLHEIGHT)
392 printf("Step: %d, proc: %d, Depth-Ave divergent at i=%d, div=%e\n",
393 prop->n,myproc,i,div_da);
394 }
395 }
396
397 mincount=0;
398 maxcount=0;
399 smin=INFTY;
400 smax=-INFTY;
401 //for(iptr=grid->celldist[0];iptr<grid->celldist[1];iptr++) {
402 for(iptr=grid->celldist[0];iptr<grid->celldist[2];iptr++) {
403 i = grid->cellp[iptr];
404
405 flag=0;
406 for(nf=0;nf<grid->nfaces[i];nf++) {
407 if(grid->mark[grid->face[i*grid->maxfaces+nf]]==2 || grid->mark[grid->face[i*grid->maxfaces+nf]]==3) {
408 flag=1;
409 break;
410 }
411 }
412
413 if(!flag) {
414 for(k=grid->ctop[i];k<grid->Nk[i];k++) {
415 if(scal[i][k]>smax) {
416 smax=scal[i][k];
417 imax=i;
418 kmax=k;
419 }
420 if(scal[i][k]<smin) {
421 smin=scal[i][k];
422 imin=i;
423 kmin=k;
424 }
425
426 if(scal[i][k]>smax_value+SMALL_CONSISTENCY && grid->dzz[i][k]>DRYCELLHEIGHT)
427 maxcount++;
428 if(scal[i][k]<smin_value-SMALL_CONSISTENCY && grid->dzz[i][k]>DRYCELLHEIGHT)
429 mincount++;
430 }
431 }
432 }
433 MPI_Reduce(&mincount,&allmincount,1,MPI_INT,MPI_SUM,0,comm);
434 MPI_Reduce(&maxcount,&allmaxcount,1,MPI_INT,MPI_SUM,0,comm);
435
436 if(mincount!=0 || maxcount!=0)
437 printf("Not CWC, step: %d, proc: %d, smin = %e at i=%d,H=%e, smax = %e at i=%d,H=%e\n",
438 prop->n,myproc,
439 smin,imin,phys->h[imin]+grid->dv[imin],
440 smax,imax,phys->h[imax]+grid->dv[imax]);
441
442 if(myproc==0 && (allmincount !=0 || allmaxcount !=0))
443 printf("Total number of CWC violations (all procs): s<s_min: %d, s>s_max: %d\n",
444 allmincount,allmaxcount);
445 }
446 }
下面介绍解读UpdateScalars函数过程：
1. 首先作为一个复杂的非静压N－S模型，变量比较多是很正常的，所以不要纠结每个变量是什么意思，能看懂就看，看不懂就猜好了。
2.要从整体入手。根据目前已知信息，这是用有限体积法求解对流扩散方程模块，而所求标量值应该就是就是第5个参数 **scal 所代表的变量。那么从函数最后一次更新 scal 变量的地方，或许能获得些许线索。
第311行：
if(grid->Nk[i]-ktop>1) TriSolve(a,b,c,d,&(scal[i][ktop]),grid->Nk[i]-ktop);
检查 TriSolve 函数的定义，原来是求解三对角方程组的解法，a，b，c 分别是系数矩阵三个对角向量，d是等号右端常向量，未知数为以 scal[i][ktop] 起始的数组。
而准备a，b，c 等系数向量时，循环变量多是按照某个三棱柱各层从上到下进行循环，所以不难看出，这个方程组求解的应该就是某个三棱柱单元体内各层标量值大小。
3. 将程序数值离散过程和理论结合起来，了解程序细节
FVCOM 模型求解过程
FVCOM 也是使用有限体积方法，但是求解和 SUNTANS 有很大不同。由于FVCOM并没有介绍对流扩散方程求解具体过程的文献，这时程序看起来就比较头疼，只能全部通过读程序来一点点理解。
FVCOM 扩散方程计算主要在程序 mod_scal.F 中，模块内全部程序如下
1 !/===========================================================================/
2 ! Copyright (c) 2007, The University of Massachusetts Dartmouth
3 ! Produced at the School of Marine Science & Technology
4 ! Marine Ecosystem Dynamics Modeling group
5 ! All rights reserved.
6 !
7 ! FVCOM has been developed by the joint UMASSD-WHOI research team. For
8 ! details of authorship and attribution of credit please see the FVCOM
9 ! technical manual or contact the MEDM group.
10 !
11 !
12 ! This file is part of FVCOM. For details, see http://fvcom.smast.umassd.edu
13 ! The full copyright notice is contained in the file COPYRIGHT located in the
14 ! root directory of the FVCOM code. This original header must be maintained
15 ! in all distributed versions.
16 !
17 ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
18 ! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
19 ! THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
20 ! PURPOSE ARE DISCLAIMED.
21 !
22 !/---------------------------------------------------------------------------/
23 ! CVS VERSION INFORMATION
24 ! $Id$
25 ! $Name$
26 ! $Revision$
27 !/===========================================================================/
28
29 !=======================================================================
30 ! FVCOM Scalar Module
31 !
32 ! contains methods:
33 ! Adv_Scal => Advect a Scalar Quantity
34 ! Vdif_Scal => Vertical Diffusion of Scalar Quantity
35 ! Bcond_Scal_OBC => Open Boundary Condition for Scalar
36 ! Bcond_Scal_PTsource => Point Sources of Scalar
37 !=======================================================================
38 Module Scalar
39
40 logical, parameter :: debug = .true.
41
42 contains
43 !==============================================================================|
44 ! Calculate Horizontal Advection and Diffusion For Scalar (f) |
45 !==============================================================================|
46 Subroutine Adv_Scal(f,fn,d_fdis,fdis,d_fflux,fflux_obc,deltat,source)
47 !------------------------------------------------------------------------------|
48
49 use all_vars
50 use lims, only: m,mt,n,nt,kbm1,kb
51 use bcs
52 use mod_obcs
53 # if defined (MULTIPROCESSOR)
54 use mod_par
55 # endif
56 # if defined (WET_DRY)
57 use mod_wd
58 # endif
59
60 # if defined (THIN_DAM)
61 use mod_dam, only : kdam,N_DAM_MATCH,IS_DAM
62 # endif
63
64 implicit none
65 real(sp), intent(in ), dimension(0:mt,kb) :: f
66 real(sp), intent(out), dimension(0:mt,kb) :: fn
67 integer , intent(in ) :: d_fdis
68 real(sp), intent(in ), dimension(d_fdis) :: fdis
69 integer , intent(in ) :: d_fflux
70 real(sp), intent(out), dimension(d_fflux,kbm1) :: fflux_obc
71 real(sp), intent(in ) :: deltat
72 logical , intent(in ) :: source
73
74 !----------------local--------------------------------------
75 real(sp), dimension(0:mt,kb) :: xflux,xflux_adv
76 real(sp), dimension(m) :: pupx,pupy,pvpx,pvpy
77 real(sp), dimension(m) :: pfpx,pfpy,pfpxd,pfpyd,viscoff
78 real(sp), dimension(3*nt) :: dtij
79 real(sp), dimension(3*nt,kbm1) :: uvn
80 real(sp), dimension(kb) :: vflux
81 real(sp) :: utmp,vtmp,sitai,ffd,ff1,x11,y11,x22,y22,x33,y33
82 real(sp) :: tmp1,tmp2,xi,yi
83 real(sp) :: dxa,dya,dxb,dyb,fij1,fij2,un
84 real(sp) :: txx,tyy,fxx,fyy,viscof,exflux,temp,fpoint
85 real(sp) :: fact,fm1,fmean
86 integer :: i,i1,i2,ia,ib,j,j1,j2,k,jtmp,jj
87 # if defined (SPHERICAL)
88 real(sp) :: ty,txpi,typi
89 # endif
90
91 # if defined (THIN_DAM)
92 INTEGER :: NX
93 real(sp) :: tmpflx
94 real(sp),dimension(kb) :: wvel
95 # endif
96
97
98 !------------------------------------------------------------------------------!
99
100 !-------------------------------------------------------
101 !Calculate Mean Values
102 !-------------------------------------------------------
103
104 fmean = sum(f(1:m,1:kbm1))/float(m*kbm1)
105
106 !-------------------------------------------------------
107 !Initialize Multipliers to Control Horizontal Diff
108 !-------------------------------------------------------
109
110 fact = 0.0_sp
111 fm1 = 1.0_sp
112 if(HORIZONTAL_MIXING_TYPE == 'closure') then
113 fact = 1.0_sp
114 fm1 = 0.0_sp
115 end if
116
117 !-------------------------------------------------------
118 !Initialize Fluxes
119 !-------------------------------------------------------
120 xflux = 0.0_sp
121 xflux_adv = 0.0_sp
122
123 !-------------------------------------------------------
124 !Calculate Normal Velocity on Control Volume Edges
125 !-------------------------------------------------------
126 !!# if !defined (WET_DRY)
127 do i=1,ncv
128 i1=ntrg(i)
129 dtij(i)=dt1(i1)
130 do k=1,kbm1
131 uvn(i,k) = v(i1,k)*dltxe(i) - u(i1,k)*dltye(i)
132
133 # if defined(PLBC)
134 uvn(i,k) = - u(i1,k)*dltye(i)
135 # endif
136
137 end do
138 end do
139 !!# else
140 !! do i=1,ncv
141 !! i1=ntrg(i)
142 !! dtij(i)=dt1(i1)
143 !! do k=1,kbm1
144 !! uvn(i,k) = vs(i1,k)*dltxe(i) - us(i1,k)*dltye(i)
145 !! end do
146 !! end do
147 !!# endif
148
149 !
150 !--Calculate the Advection and Horizontal Diffusion Terms----------------------!
151 !
152
153 do k=1,kbm1
154 pfpx = 0.0_sp
155 pfpy = 0.0_sp
156 pfpxd = 0.0_sp
157 pfpyd = 0.0_sp
158 do i=1,m
159 do j=1,ntsn(i)-1
160 i1=nbsn(i,j)
161 i2=nbsn(i,j+1)
162
163 # if defined (WET_DRY)
164 IF(ISWETN(I1) == 0 .AND. ISWETN(I2) == 1)THEN
165 FFD=0.5_SP*(f(I,K)+f(I2,K))
166 FF1=0.5_SP*(f(I,K)+f(I2,K))
167 ELSE IF(ISWETN(I1) == 1 .AND. ISWETN(I2) == 0)THEN
168 FFD=0.5_SP*(f(I1,K)+f(I,K))
169 FF1=0.5_SP*(f(I1,K)+f(I,K))
170 ELSE IF(ISWETN(I1) == 0 .AND. ISWETN(I2) == 0)THEN
171 FFD=0.5_SP*(f(I,K)+f(I,K))
172 FF1=0.5_SP*(f(I,K)+f(I,K))
173 ELSE
174 FFD=0.5_SP*(f(I1,K)+f(I2,K))
175 FF1=0.5_SP*(f(I1,K)+f(I2,K))
176 END IF
177 # else
178 ffd=0.5_sp*(f(i1,k)+f(i2,k)) !-fmean1(i1,k)-fmean1(i2,k))
179 ff1=0.5_sp*(f(i1,k)+f(i2,k))
180 # endif
181
182 # if defined (SPHERICAL)
183 ty=0.5_sp*(vy(i1)+vy(i2))
184 txpi=(vx(i2)-vx(i1))*tpi*cos(deg2rad*ty)
185 typi=(vy(i1)-vy(i2))*tpi
186 pfpx(i)=pfpx(i)+ff1*typi
187 pfpy(i)=pfpy(i)+ff1*txpi
188 pfpxd(i)=pfpxd(i)+ffd*typi
189 pfpyd(i)=pfpyd(i)+ffd*txpi
190 # else
191 pfpx(i) = pfpx(i) +ff1*(vy(i1)-vy(i2))
192 pfpy(i) = pfpy(i) +ff1*(vx(i2)-vx(i1))
193 pfpxd(i)= pfpxd(i)+ffd*(vy(i1)-vy(i2))
194 pfpyd(i)= pfpyd(i)+ffd*(vx(i2)-vx(i1))
195 # endif
196 end do
197
198 ! gather all neighboring control volumes connecting at dam node
199 # if defined (THIN_DAM)
200 IF(IS_DAM(I)==1.AND.K<=KDAM(I))THEN
201 DO NX=1,N_DAM_MATCH(I,1)
202 DO J=1,NTSN(N_DAM_MATCH(I,NX+1))-1
203 I1=NBSN(N_DAM_MATCH(I,NX+1),J)
204 I2=NBSN(N_DAM_MATCH(I,NX+1),J+1)
205 FFD=0.5_SP*(f(I1,K)+f(I2,K)) !-SMEAN1(I1,K)-SMEAN1(I2,K))
206 FF1=0.5_SP*(f(I1,K)+f(I2,K))
207 # if defined (SPHERICAL)
208 ty=0.5_sp*(vy(i1)+vy(i2))
209 txpi=(vx(i2)-vx(i1))*tpi*cos(deg2rad*ty)
210 typi=(vy(i1)-vy(i2))*tpi
211 pfpx(i)=pfpx(i)+ff1*typi
212 pfpy(i)=pfpy(i)+ff1*txpi
213 pfpxd(i)=pfpxd(i)+ffd*typi
214 pfpyd(i)=pfpyd(i)+ffd*txpi
215 # else
216 pfpx(i) = pfpx(i) +ff1*(vy(i1)-vy(i2))
217 pfpy(i) = pfpy(i) +ff1*(vx(i2)-vx(i1))
218 pfpxd(i)= pfpxd(i)+ffd*(vy(i1)-vy(i2))
219 pfpyd(i)= pfpyd(i)+ffd*(vx(i2)-vx(i1))
220 # endif
221 END DO
222 END DO
223 END IF
224 # endif
225
226 # if !defined (THIN_DAM)
227 pfpx(i) =pfpx(i )/art2(i)
228 pfpy(i) =pfpy(i )/art2(i)
229 pfpxd(i) =pfpxd(i)/art2(i)
230 pfpyd(i) =pfpyd(i)/art2(i)
231 # else
232 IF(IS_DAM(I)==1.AND.K<=KDAM(I))THEN
233 PFPX(I)=PFPX(I)/(ART2(I)+SUM(ART2(N_DAM_MATCH(I,2:1+N_DAM_MATCH(I,1)))))
234 PFPY(I)=PFPY(I)/(ART2(I)+SUM(ART2(N_DAM_MATCH(I,2:1+N_DAM_MATCH(I,1)))))
235 PFPXD(I)=PFPXD(I)/(ART2(I)+SUM(ART2(N_DAM_MATCH(I,2:1+N_DAM_MATCH(I,1)))))
236 PFPYD(I)=PFPYD(I)/(ART2(I)+SUM(ART2(N_DAM_MATCH(I,2:1+N_DAM_MATCH(I,1)))))
237 ELSE
238 PFPX(I)=PFPX(I)/ART2(I)
239 PFPY(I)=PFPY(I)/ART2(I)
240 PFPXD(I)=PFPXD(I)/ART2(I)
241 PFPYD(I)=PFPYD(I)/ART2(I)
242 END IF
243 # endif
244
245 end do
246
247 if(k == kbm1)then
248 do i=1,m
249 pfpxb(i) = pfpx(i)
250 pfpyb(i) = pfpy(i)
251 end do
252 end if
253
254 do i=1,m
255 pupx(i)=0.0_sp
256 pupy(i)=0.0_sp
257 pvpx(i)=0.0_sp
258 pvpy(i)=0.0_sp
259 j=1
260 i1=nbve(i,j)
261 jtmp=nbvt(i,j)
262 j1=jtmp+1-(jtmp+1)/4*3
263 j2=jtmp+2-(jtmp+2)/4*3
264 x11=0.5_sp*(vx(i)+vx(nv(i1,j1)))
265 y11=0.5_sp*(vy(i)+vy(nv(i1,j1)))
266 x22=xc(i1)
267 y22=yc(i1)
268 x33=0.5_sp*(vx(i)+vx(nv(i1,j2)))
269 y33=0.5_sp*(vy(i)+vy(nv(i1,j2)))
270
271 # if defined (SPHERICAL)
272 ty =0.5_sp*(y11+y33)
273 txpi=(x33-x11)*tpi*cos(deg2rad*ty)
274 typi=(y11-y33)*tpi
275 pupx(i)=pupx(i)+u(i1,k)*typi
276 pupy(i)=pupy(i)+u(i1,k)*txpi
277 pvpx(i)=pvpx(i)+v(i1,k)*typi
278 pvpy(i)=pvpy(i)+v(i1,k)*txpi
279 # else
280 pupx(i)=pupx(i)+u(i1,k)*(y11-y33)
281 pupy(i)=pupy(i)+u(i1,k)*(x33-x11)
282 pvpx(i)=pvpx(i)+v(i1,k)*(y11-y33)
283 pvpy(i)=pvpy(i)+v(i1,k)*(x33-x11)
284 # endif
285
286 if(isonb(i) /= 0) then
287 # if defined (SPHERICAL)
288 ty=0.5_sp*(vy(i)+y11)
289 txpi=(x11-vx(i))*tpi*cos(deg2rad*ty)
290 typi=(vy(i)-y11)*tpi
291 pupx(i)=pupx(i)+u(i1,k)*typi
292 pupy(i)=pupy(i)+u(i1,k)*txpi
293 pvpx(i)=pvpx(i)+v(i1,k)*typi
294 pvpy(i)=pvpy(i)+v(i1,k)*txpi
295 # else
296 pupx(i)=pupx(i)+u(i1,k)*(vy(i)-y11)
297 pupy(i)=pupy(i)+u(i1,k)*(x11-vx(i))
298 pvpx(i)=pvpx(i)+v(i1,k)*(vy(i)-y11)
299 pvpy(i)=pvpy(i)+v(i1,k)*(x11-vx(i))
300 # endif
301 end if
302
303 do j=2,ntve(i)-1
304 i1=nbve(i,j)
305 jtmp=nbvt(i,j)
306 j1=jtmp+1-(jtmp+1)/4*3
307 j2=jtmp+2-(jtmp+2)/4*3
308 x11=0.5_sp*(vx(i)+vx(nv(i1,j1)))
309 y11=0.5_sp*(vy(i)+vy(nv(i1,j1)))
310 x22=xc(i1)
311 y22=yc(i1)
312 x33=0.5_sp*(vx(i)+vx(nv(i1,j2)))
313 y33=0.5_sp*(vy(i)+vy(nv(i1,j2)))
314
315 # if defined (SPHERICAL)
316 ty=0.5_sp*(y11+y33)
317 txpi=(x33-x11)*tpi*COS(deg2rad*TY)
318 typi=(y11-y33)*tpi
319 pupx(i)=pupx(i)+u(i1,k)*typi
320 pupy(i)=pupy(i)+u(i1,k)*txpi
321 pvpx(i)=pvpx(i)+v(i1,k)*typi
322 pvpy(i)=pvpy(i)+v(i1,k)*txpi
323 # else
324 pupx(i)=pupx(i)+u(i1,k)*(y11-y33)
325 pupy(i)=pupy(i)+u(i1,k)*(x33-x11)
326 pvpx(i)=pvpx(i)+v(i1,k)*(y11-y33)
327 pvpy(i)=pvpy(i)+v(i1,k)*(x33-x11)
328 # endif
329 end do
330 j=ntve(i)
331 i1=nbve(i,j)
332 jtmp=nbvt(i,j)
333 j1=jtmp+1-(jtmp+1)/4*3
334 j2=jtmp+2-(jtmp+2)/4*3
335 x11=0.5_sp*(vx(i)+vx(nv(i1,j1)))
336 y11=0.5_sp*(vy(i)+vy(nv(i1,j1)))
337 x22=xc(i1)
338 y22=yc(i1)
339 x33=0.5_sp*(vx(i)+vx(nv(i1,j2)))
340 y33=0.5_sp*(vy(i)+vy(nv(i1,j2)))
341
342 # if defined (SPHERICAL)
343 ty=0.5*(Y11+Y33)
344 txpi=(x33-x11)*tpi*cos(deg2rad*TY)
345 typi=(y11-y33)*tpi
346 pupx(i)=pupx(i)+u(i1,k)*typi
347 pupy(i)=pupy(i)+u(i1,k)*txpi
348 pvpx(i)=pvpx(i)+v(i1,k)*typi
349 pvpy(i)=pvpy(i)+v(i1,k)*txpi
350 # else
351 pupx(i)=pupx(i)+u(i1,k)*(y11-y33)
352 pupy(i)=pupy(i)+u(i1,k)*(x33-x11)
353 pvpx(i)=pvpx(i)+v(i1,k)*(y11-y33)
354 pvpy(i)=pvpy(i)+v(i1,k)*(x33-x11)
355 # endif
356
357 if(isonb(i) /= 0) then
358 # if defined (SPHERICAL)
359 ty=0.5*(Y11+VY(I))
360 txpi=(VX(I)-X11)*tpi*COS(deg2rad*ty)
361 typi=(Y11-VY(I))*tpi
362 pupx(i)=pupx(i)+u(i1,k)*typi
363 pupy(i)=pupy(i)+u(i1,k)*txpi
364 pvpx(i)=pvpx(i)+v(i1,k)*typi
365 pvpy(i)=pvpy(i)+v(i1,k)*txpi
366 # else
367 pupx(i)=pupx(i)+u(i1,k)*(y11-vy(i))
368 pupy(i)=pupy(i)+u(i1,k)*(vx(i)-x11)
369 pvpx(i)=pvpx(i)+v(i1,k)*(y11-vy(i))
370 pvpy(i)=pvpy(i)+v(i1,k)*(vx(i)-x11)
371 # endif
372 end if
373 pupx(i)=pupx(i)/art1(i)
374 pupy(i)=pupy(i)/art1(i)
375 pvpx(i)=pvpx(i)/art1(i)
376 pvpy(i)=pvpy(i)/art1(i)
377 tmp1=pupx(i)**2+pvpy(i)**2
378 tmp2=0.5_sp*(pupy(i)+pvpx(i))**2
379 viscoff(i)=sqrt(tmp1+tmp2)*art1(i)
380 end do
381 ! if(k == kbm1) then
382 ! ah_bottom(1:m) = horcon*(fact*viscoff(1:m) + fm1)
383 ! end if
384
385
386 do i=1,ncv_i
387 ia=niec(i,1)
388 ib=niec(i,2)
389 xi=0.5_sp*(xije(i,1)+xije(i,2))
390 yi=0.5_sp*(yije(i,1)+yije(i,2))
391 # if defined (SPHERICAL)
392 ty=0.5_sp*(yi+vy(ia))
393 dxa=(xi-vx(ia))*tpi*cos(deg2rad*ty)
394 dya=(yi-vy(ia))*tpi
395 ty=0.5*(YI+VY(IB))
396 DXB=(XI-VX(IB))*tpi*COS(deg2rad*ty)
397 DYB=(YI-VY(IB))*tpi
398 # else
399 dxa=xi-vx(ia)
400 dya=yi-vy(ia)
401 dxb=xi-vx(ib)
402 dyb=yi-vy(ib)
403 # endif
404 fij1=f(ia,k)+dxa*pfpx(ia)+dya*pfpy(ia)
405 fij2=f(ib,k)+dxb*pfpx(ib)+dyb*pfpy(ib)
406 un=uvn(i,k)
407
408 ! viscof=horcon*(fact*(viscoff(ia)+viscoff(ib))*0.5_sp + fm1)
409 VISCOF=(FACT*0.5_SP*(VISCOFF(IA)*NN_HVC(IA)+VISCOFF(IB)*NN_HVC(IB)) + FM1*0.5_SP*(NN_HVC(IA)+NN_HVC(IB)))
410
411 txx=0.5_sp*(pfpxd(ia)+pfpxd(ib))*viscof
412 tyy=0.5_sp*(pfpyd(ia)+pfpyd(ib))*viscof
413
414 fxx=-dtij(i)*txx*dltye(i)
415 fyy= dtij(i)*tyy*dltxe(i)
416
417 # if defined (PLBC)
418 fyy=0.0_SP
419 # endif
420
421 exflux=-un*dtij(i)* &
422 ((1.0_sp+sign(1.0_sp,un))*fij2+(1.0_sp-sign(1.0_sp,un))*fij1)*0.5_sp+fxx+fyy
423
424 xflux(ia,k)=xflux(ia,k)+exflux
425 xflux(ib,k)=xflux(ib,k)-exflux
426
427 xflux_adv(ia,k)=xflux_adv(ia,k)+(exflux-fxx-fyy)
428 xflux_adv(ib,k)=xflux_adv(ib,k)-(exflux-fxx-fyy)
429
430 # if defined (THIN_DAM)
431 IF(K<=KDAM(IA).AND.IS_DAM(IA)==1)THEN
432 IF(N_DAM_MATCH(IA,1)==1)THEN
433 XFLUX(N_DAM_MATCH(IA,2),K) = XFLUX(N_DAM_MATCH(IA,2),K) + EXFLUX
434 XFLUX_ADV(N_DAM_MATCH(IA,2),K) = XFLUX_ADV(N_DAM_MATCH(IA,2),K) +(EXFLUX-FXX-FYY)
435 END IF
436 IF(N_DAM_MATCH(IA,1)==2)THEN
437 XFLUX(N_DAM_MATCH(IA,2),K) = XFLUX(N_DAM_MATCH(IA,2),K) + EXFLUX
438 XFLUX(N_DAM_MATCH(IA,3),K) = XFLUX(N_DAM_MATCH(IA,3),K) + EXFLUX
439 XFLUX_ADV(N_DAM_MATCH(IA,2),K) = XFLUX_ADV(N_DAM_MATCH(IA,2),K) +(EXFLUX-FXX-FYY)
440 XFLUX_ADV(N_DAM_MATCH(IA,3),K) = XFLUX_ADV(N_DAM_MATCH(IA,3),K) +(EXFLUX-FXX-FYY)
441 END IF
442 IF(N_DAM_MATCH(IA,1)==3)THEN
443 XFLUX(N_DAM_MATCH(IA,2),K) = XFLUX(N_DAM_MATCH(IA,2),K) + EXFLUX
444 XFLUX(N_DAM_MATCH(IA,3),K) = XFLUX(N_DAM_MATCH(IA,3),K) + EXFLUX
445 XFLUX(N_DAM_MATCH(IA,4),K) = XFLUX(N_DAM_MATCH(IA,4),K) + EXFLUX
446 XFLUX_ADV(N_DAM_MATCH(IA,2),K) = XFLUX_ADV(N_DAM_MATCH(IA,2),K) +(EXFLUX-FXX-FYY)
447 XFLUX_ADV(N_DAM_MATCH(IA,3),K) = XFLUX_ADV(N_DAM_MATCH(IA,3),K) +(EXFLUX-FXX-FYY)
448 XFLUX_ADV(N_DAM_MATCH(IA,4),K) = XFLUX_ADV(N_DAM_MATCH(IA,4),K) +(EXFLUX-FXX-FYY)
449 END IF
450 END IF
451 IF(K<=KDAM(IB).AND.IS_DAM(IB)==1)THEN
452 IF(N_DAM_MATCH(IB,1)==1)THEN
453 XFLUX(N_DAM_MATCH(IB,2),K) = XFLUX(N_DAM_MATCH(IB,2),K) - EXFLUX
454 XFLUX_ADV(N_DAM_MATCH(IB,2),K) = XFLUX_ADV(N_DAM_MATCH(IB,2),K) - (EXFLUX-FXX-FYY)
455 END IF
456 IF(N_DAM_MATCH(IB,1)==2)THEN
457 XFLUX(N_DAM_MATCH(IB,2),K) = XFLUX(N_DAM_MATCH(IB,2),K) - EXFLUX
458 XFLUX(N_DAM_MATCH(IB,3),K) = XFLUX(N_DAM_MATCH(IB,3),K) - EXFLUX
459 XFLUX_ADV(N_DAM_MATCH(IB,2),K) = XFLUX_ADV(N_DAM_MATCH(IB,2),K) - (EXFLUX-FXX-FYY)
460 XFLUX_ADV(N_DAM_MATCH(IB,3),K) = XFLUX_ADV(N_DAM_MATCH(IB,3),K) - (EXFLUX-FXX-FYY)
461 END IF
462 IF(N_DAM_MATCH(IB,1)==3)THEN
463 XFLUX(N_DAM_MATCH(IB,2),K) = XFLUX(N_DAM_MATCH(IB,2),K) - EXFLUX
464 XFLUX(N_DAM_MATCH(IB,3),K) = XFLUX(N_DAM_MATCH(IB,3),K) - EXFLUX
465 XFLUX(N_DAM_MATCH(IB,4),K) = XFLUX(N_DAM_MATCH(IB,4),K) - EXFLUX
466 XFLUX_ADV(N_DAM_MATCH(IB,2),K) = XFLUX_ADV(N_DAM_MATCH(IB,2),K) - (EXFLUX-FXX-FYY)
467 XFLUX_ADV(N_DAM_MATCH(IB,3),K) = XFLUX_ADV(N_DAM_MATCH(IB,3),K) - (EXFLUX-FXX-FYY)
468 XFLUX_ADV(N_DAM_MATCH(IB,4),K) = XFLUX_ADV(N_DAM_MATCH(IB,4),K) - (EXFLUX-FXX-FYY)
469 END IF
470 END IF
471 # endif
472 end do
473 end do !!sigma loop
474
475 !---------------------------------------------------------------------------------
476 ! Accumulate Fluxes at Boundary Nodes
477 !---------------------------------------------------------------------------------
478
479 # if defined (MULTIPROCESSOR)
480 if(par)call node_match(0,nbn,bn_mlt,bn_loc,bnc,mt,kb,myid,nprocs,xflux,xflux_adv)
481 # endif
482
483 !---------------------------------------------------------------------------------
484 ! Store Advective Fluxes at the Boundary
485 !---------------------------------------------------------------------------------
486 do k=1,kbm1
487 if(iobcn > 0) then
488 do i=1,iobcn
489 i1=i_obc_n(i)
490 fflux_obc(i,k)=xflux_adv(i1,k)
491 end do
492 end if
493 end do
494
495 !---------------------------------------------------------------------------------
496 ! Calculate Vertical Advection Terms
497 !---------------------------------------------------------------------------------
498
499 do i=1,m
500 # if defined (WET_DRY)
501 if(iswetn(i)*iswetnt(i) == 1) then
502 # endif
503 # if defined (THIN_DAM)
504 if(IS_DAM(I)==1)then
505 wvel(1:kb)=0.0_sp
506 call calc_vflux(kbm1,f(i,1:kbm1),wvel(1:kb),vflux)
507 else
508 call calc_vflux(kbm1,f(i,1:kbm1),wts(i,1:kb),vflux)
509 end if
510 # else
511 call calc_vflux(kbm1,f(i,1:kbm1),wts(i,1:kb),vflux)
512 # endif
513
514 do k=1,kbm1
515 if(isonb(i) == 2) then
516 xflux(i,k)= (vflux(k)-vflux(k+1))*art1(i)/dz(i,k)
517 else
518 xflux(i,k)=xflux(i,k)+ (vflux(k)-vflux(k+1))*art1(i)/dz(i,k)
519 end if
520 # if defined (THIN_DAM)
521 IF(IS_DAM(I)==1.AND.K<=KDAM(I))THEN
522 tmpflx = (vflux(k)-vflux(k+1))*art1(i)/dz(i,k)
523 IF(N_DAM_MATCH(I,1)==1)THEN
524 XFLUX(N_DAM_MATCH(I,2),K) = XFLUX(N_DAM_MATCH(I,2),K)+tmpflx
525 END IF
526 IF(N_DAM_MATCH(I,1)==2)THEN
527 XFLUX(N_DAM_MATCH(I,2),K) = XFLUX(N_DAM_MATCH(I,2),K)+tmpflx
528 XFLUX(N_DAM_MATCH(I,3),K) = XFLUX(N_DAM_MATCH(I,3),K)+tmpflx
529 END IF
530 IF(N_DAM_MATCH(I,1)==3)THEN
531 XFLUX(N_DAM_MATCH(I,2),K) = XFLUX(N_DAM_MATCH(I,2),K)+tmpflx
532 XFLUX(N_DAM_MATCH(I,3),K) = XFLUX(N_DAM_MATCH(I,3),K)+tmpflx
533 XFLUX(N_DAM_MATCH(I,4),K) = XFLUX(N_DAM_MATCH(I,4),K)+tmpflx
534 END IF
535 END IF
536 # endif
537 end do
538 # if defined (WET_DRY)
539 end if
540 # endif
541 end do
542
543 !-------------------------------------------------------
544 !Point Source
545 !-------------------------------------------------------
546 if(source)then !!user specified
547
548 if(RIVER_TS_SETTING == 'calculated') then
549 if(RIVER_INFLOW_LOCATION == 'node') then
550 do j=1,numqbc
551 jj=inodeq(j)
552 fpoint=fdis(j)
553 do k=1,kbm1
554 xflux(jj,k)=xflux(jj,k) - qdis(j)*vqdist(j,k)*fpoint !/dz(jj,k)
555 end do
556 end do
557 else if(RIVER_INFLOW_LOCATION == 'edge') then
558 write(*,*)'scalar advection not setup for "edge" point source'
559 stop
560 end if
561 end if
562
563 else
564
565 if(RIVER_TS_SETTING == 'calculated') then
566 if(RIVER_INFLOW_LOCATION == 'node') then
567 do j=1,numqbc
568 jj=inodeq(j)
569 do k=1,kbm1
570 fpoint = f(jj,k)
571 xflux(jj,k)=xflux(jj,k) - qdis(j)*vqdist(j,k)*fpoint !/dz(jj,k)
572 end do
573 end do
574 else if(RIVER_INFLOW_LOCATION == 'edge') then
575 write(*,*)'scalar advection not setup for "edge" point source'
576 stop
577 end if
578 end if
579
580 endif
581 !------------------------------------------------------------------------
582 !Update Scalar Quantity
583 !------------------------------------------------------------------------
584
585 do i=1,m
586 # if defined (WET_DRY)
587 if(iswetn(i)*iswetnt(i) == 1 )then
588 # endif
589 do k=1,kbm1
590 # if !defined (THIN_DAM)
591 fn(i,k)=(f(i,k)-xflux(i,k)/art1(i)*(deltat/dt(i)))*(dt(i)/dtfa(i))
592 # else
593 IF(IS_DAM(I)==1.AND.K<=KDAM(I))THEN
594 fn(i,k)=(f(i,k)-xflux(i,k)/(ART1(I)&
595 &+SUM(ART1(N_DAM_MATCH(I,2:1+N_DAM_MATCH(I,1)))))*(deltat/dt(i)))*(dt(i)/dtfa(i))
596 ELSE
597 fn(i,k)=(f(i,k)-xflux(i,k)/art1(i)*(deltat/dt(i)))*(dt(i)/dtfa(i))
598 END IF
599 # endif
600 end do
601 # if defined (WET_DRY)
602 else
603 do k=1,kbm1
604 fn(i,k)=f(i,k)
605 end do
606 end if
607 # endif
608 end do
609
610 return
611 End Subroutine Adv_Scal
612 !==============================================================================|
613
614 !==============================================================================|
615 ! Vertical Diffusion of Scalar |
616 !==============================================================================|
617 Subroutine Vdif_Scal(f,deltat)
618
619 use mTridiagonal
620 use all_vars
621 # if defined (THIN_DAM)
622 use mod_dam,only : NODE_DAM1_N,NODE_DAM2_N,NODE_DAM3_N, &
623 &I_NODE_DAM1_N,I_NODE_DAM2_N,I_NODE_DAM3_N, &
624 &kdam
625 # endif
626
627 Implicit None
628 Real(sp), intent(inout) :: f(0:mt,kb)
629 Real(sp), intent(in ) :: deltat
630 !--local--------------------
631 integer :: i,k,ll
632 real(sp) :: dsqrd,dfdz,visb
633 real(sp) :: fsol(0:kb)
634
635 # if defined (THIN_DAM)
636 real(sp) :: ftmp,stmp
637 # endif
638
639 call init_tridiagonal(kb)
640
641 Do i=1,m
642 dsqrd = d(i)*d(i)
643
644 !----------------------------------------------------------------
645 ! Set up Diagonals of Matrix (lower=au,diag=bu,upper=cu)
646 !----------------------------------------------------------------
647
648
649 !Surface
650 au(1) = 0.0
651 cu(1)= - deltat*(kh(i,2)+umol)/(dzz(i,1)*dz(i,1)*dsqrd)
652 bu(1)= 1.0 - cu(1)
653
654 !Interior
655 do k=2,kbm1-1
656 au(k) = - deltat*(kh(i,k )+umol)/(dzz(i,k-1)*dz(i,k)*dsqrd)
657 cu(k) = - deltat*(kh(i,k+1)+umol)/(dzz(i,k )*dz(i,k)*dsqrd)
658 bu(k) = 1.0 - cu(k) - au(k)
659 end do
660
661 !Bottom
662 au(kbm1) = - deltat*(kh(i,kbm1)+umol)/(dzz(i,kbm1-1)*dz(i,kbm1)*dsqrd)
663 cu(kbm1) = 0.0
664 bu(kbm1) = 1.0 - au(kbm1)
665
666 !----------------------------------------------------------------
667 ! Set up RHS forcing vector and boundary conditions
668 !----------------------------------------------------------------
669 do k=1,kbm1
670 du(k) = f(i,k)
671 end do
672
673 !Free Surface: No flux
674
675 !Bottom: No flux
676
677
678 !----------------------------------------------------------------
679 ! Solve
680 !----------------------------------------------------------------
681
682 call tridiagonal(kb,1,kbm1,fsol)
683
684 !Transfer
685 f(i,1:kbm1) = fsol(1:kbm1)
686
687 End Do
688
689 # if defined (THIN_DAM)
690 DO K=1,KBM1
691 DO I=1,NODE_DAM1_N
692 IF(K<=KDAM(I_NODE_DAM1_N(I,1)).AND.K<=KDAM(I_NODE_DAM1_N(I,2)) )THEN
693 FTMP=F(I_NODE_DAM1_N(I,1),K)*ART1(I_NODE_DAM1_N(I,1)) &
694 & +F(I_NODE_DAM1_N(I,2),K)*ART1(I_NODE_DAM1_N(I,2))
695 STMP=ART1(I_NODE_DAM1_N(I,1))+ART1(I_NODE_DAM1_N(I,2))
696 F(I_NODE_DAM1_N(I,1),K)=FTMP/STMP
697 F(I_NODE_DAM1_N(I,2),K)=FTMP/STMP
698 END IF
699 END DO
700
701 DO I=1,NODE_DAM2_N
702 IF(K<=KDAM(I_NODE_DAM2_N(I,1)).AND.K<=KDAM(I_NODE_DAM2_N(I,2)) &
703 & .AND.K<=KDAM(I_NODE_DAM2_N(I,2)) )THEN
704 FTMP= F(I_NODE_DAM2_N(I,1),K)*ART1(I_NODE_DAM2_N(I,1)) &
705 & +F(I_NODE_DAM2_N(I,2),K)*ART1(I_NODE_DAM2_N(I,2)) &
706 & +F(I_NODE_DAM2_N(I,3),K)*ART1(I_NODE_DAM2_N(I,3))
707 STMP=ART1(I_NODE_DAM2_N(I,1))+ART1(I_NODE_DAM2_N(I,2)) &
708 & +ART1(I_NODE_DAM2_N(I,3))
709 F(I_NODE_DAM2_N(I,1),K)=FTMP/STMP
710 F(I_NODE_DAM2_N(I,2),K)=FTMP/STMP
711 F(I_NODE_DAM2_N(I,3),K)=FTMP/STMP
712 END IF
713 END DO
714
715 DO I=1,NODE_DAM3_N
716 IF(K<=KDAM(I_NODE_DAM3_N(I,1)).AND.K<=KDAM(I_NODE_DAM3_N(I,2)) &
717 & .AND.K<=KDAM(I_NODE_DAM3_N(I,3)).AND.K<=KDAM(I_NODE_DAM3_N(I,4)) )THEN
718 FTMP =F(I_NODE_DAM3_N(I,1),K)*ART1(I_NODE_DAM3_N(I,1)) &
719 & +F(I_NODE_DAM3_N(I,2),K)*ART1(I_NODE_DAM3_N(I,2)) &
720 & +F(I_NODE_DAM3_N(I,3),K)*ART1(I_NODE_DAM3_N(I,3)) &
721 & +F(I_NODE_DAM3_N(I,4),K)*ART1(I_NODE_DAM3_N(I,4))
722 STMP =ART1(I_NODE_DAM3_N(I,1)) + ART1(I_NODE_DAM3_N(I,2)) &
723 & + ART1(I_NODE_DAM3_N(I,3)) + ART1(I_NODE_DAM3_N(I,4))
724 F(I_NODE_DAM3_N(I,1),K)=FTMP/STMP
725 F(I_NODE_DAM3_N(I,2),K)=FTMP/STMP
726 F(I_NODE_DAM3_N(I,3),K)=FTMP/STMP
727 F(I_NODE_DAM3_N(I,4),K)=FTMP/STMP
728 END IF
729 END DO
730 END DO
731 # endif
732
733
734 End Subroutine Vdif_Scal
735
736
737 !==============================================================================|
738 ! Set Point Source Conditions for Scalar Function |
739 !==============================================================================|
740
741 Subroutine Bcond_Scal_PTsource(f,fn,fdis)
742
743 !------------------------------------------------------------------------------|
744 use all_vars
745 use bcs
746 use mod_obcs
747 implicit none
748 real(sp), intent(in ), dimension(0:mt,kb) :: f
749 real(sp), intent(out), dimension(0:mt,kb) :: fn
750 real(sp), intent(in ), dimension(numqbc ) :: fdis
751 !--local-------------------------------------------
752 integer :: i,j,k,j1,j11,j22
753 !------------------------------------------------------------------------------|
754
755
756 !--------------------------------------------
757 ! Set Source Terms
758 !--------------------------------------------
759 if(RIVER_TS_SETTING == 'specified') then
760 if(numqbc > 0) then
761 if(RIVER_INFLOW_LOCATION == 'node') then
762 do i=1,numqbc
763 j11=inodeq(i)
764 do k=1,kbm1
765 fn(j11,k)=fdis(i)
766 end do
767 end do
768 else if(RIVER_INFLOW_LOCATION == 'edge') then
769 do i=1,numqbc
770 j11=n_icellq(i,1)
771 j22=n_icellq(i,2)
772 do k=1,kbm1
773 fn(j11,k)=fdis(i)
774 fn(j22,k)=fdis(i)
775 end do
776 end do
777 end if
778 end if
779 end if
780
781 return
782 End Subroutine Bcond_Scal_PTSource
783 !==============================================================================|
784 !==============================================================================|
785
786 !==============================================================================|
787 ! Set Boundary Conditions for Scalar Function on Open Boundary |
788 !==============================================================================|
789
790 Subroutine Bcond_Scal_OBC(f,fn,fflux_obc,f_obc,deltat,alpha_nudge)
791
792 !------------------------------------------------------------------------------|
793 use all_vars
794 use bcs
795 use mod_obcs
796 implicit none
797 real(sp), intent(in ), dimension(0:mt,kb) :: f
798 real(sp), intent(inout), dimension(0:mt,kb) :: fn
799 real(sp), intent(in ), dimension(iobcn+1,kbm1) :: fflux_obc
800 real(sp), intent(in ), dimension(iobcn ) :: f_obc
801 real(sp), intent(in ) :: deltat
802 real(sp), intent(in ) :: alpha_nudge
803 !--local-------------------------------------------
804 real(sp) :: f2d,f2d_next,f2d_obc,xflux2d,tmp
805 integer :: i,j,k,j1,j11,j22
806 !------------------------------------------------------------------------------|
807
808 !--------------------------------------------
809 ! Set Scalar Value on Open Boundary
810 !--------------------------------------------
811 if(iobcn > 0) then
812 do i=1,iobcn
813 j=i_obc_n(i)
814 j1=next_obc(i)
815 f2d=0.0_sp
816 f2d_next=0.0_sp
817 xflux2d=0.0_sp
818 do k=1,kbm1
819 f2d=f2d+f(j,k)*dz(j,k)
820 f2d_next=f2d_next+fn(j1,k)*dz(j1,k)
821 xflux2d=xflux2d+fflux_obc(i,k)*dz(j,k)
822 end do
823
824 if(uard_obcn(i) > 0.0_sp) then
825 tmp=xflux2d+f2d*uard_obcn(i)
826 f2d_obc=(f2d*dt(j)-tmp*deltat/art1(j))/d(j)
827 do k=1,kbm1
828 fn(j,k)=fn(j1,k) !f2d_obc+(fn(j1,k)-f2d_next)
829 end do
830 else
831 do k=1,kbm1
832 fn(j,k) = f(j,k)-alpha_nudge*(f(j,k)-f_obc(i))
833 end do
834 end if
835 end do
836 endif
837
838 return
839 End Subroutine Bcond_Scal_OBC
840 !==============================================================================|
841 !==============================================================================|
842
843 Subroutine fct_sed(f,fn)
844 !==============================================================================|
845 USE ALL_VARS
846 USE MOD_UTILS
847 USE BCS
848 USE MOD_OBCS
849 IMPLICIT NONE
850 real(sp), intent(inout), dimension(0:mt,kb) :: fn
851 real(sp), intent(in), dimension(0:mt,kb) :: f
852 REAL(SP):: SMAX,SMIN
853 INTEGER :: I,J,K,K1
854 !==============================================================================|
855 IF(DBG_SET(DBG_SBR)) WRITE(IPT,*)"Start: fct_sed"
856
857 nodes: DO I=1,M
858
859 ! SKIP OPEN BOUNDARY NODES
860 IF(IOBCN > 0)THEN
861 DO J=1,IOBCN
862 IF(I == I_OBC_N(J)) CYCLE nodes
863 END DO
864 END IF
865
866 ! SKIP RIVER INFLOW POINTS
867 IF(NUMQBC > 0)THEN
868 DO J=1,NUMQBC
869 IF(RIVER_INFLOW_LOCATION == 'node')THEN
870 IF(I == INODEQ(J)) CYCLE nodes
871 END IF
872 IF(RIVER_INFLOW_LOCATION == 'edge')THEN
873 IF(I == N_ICELLQ(J,1) .OR. I == N_ICELLQ(J,2)) CYCLE nodes
874 END IF
875 END DO
876 END IF
877
878 ! SKIP GROUND WATER INFLOW POINTS
879 IF(BFWDIS(I) .GT. 0.0_SP .and. GROUNDWATER_SALT_ON) CYCLE nodes
880
881 K1 = 1
882 IF(PRECIPITATION_ON) K1 = 2
883 ! DO K=1,KBM1
884 DO K=K1,KBM1
885 SMAX = MAXVAL(f(NBSN(I,1:NTSN(I)),K))
886 SMIN = MINVAL(f(NBSN(I,1:NTSN(I)),K))
887
888 IF(K == 1)THEN
889 SMAX = MAX(SMAX,(f(I,K)*DZ(I,K+1)+f(I,K+1)*DZ(I,K))/ &
890 (DZ(I,K)+DZ(I,K+1)))
891 SMIN = MIN(SMIN,(f(I,K)*DZ(I,K+1)+f(I,K+1)*DZ(I,K))/ &
892 (DZ(I,K)+DZ(I,K+1)))
893 ELSE IF(K == KBM1)THEN
894 SMAX = MAX(SMAX,(f(I,K)*DZ(I,K-1)+f(I,K-1)*DZ(I,K))/ &
895 (DZ(I,K)+DZ(I,K-1)))
896 SMIN = MIN(SMIN,(f(I,K)*DZ(I,K-1)+f(I,K-1)*DZ(I,K))/ &
897 (DZ(I,K)+DZ(I,K-1)))
898 ELSE
899 SMAX = MAX(SMAX,(f(I,K)*DZ(I,K-1)+f(I,K-1)*DZ(I,K))/ &
900 (DZ(I,K)+DZ(I,K-1)), &
901 (f(I,K)*DZ(I,K+1)+f(I,K+1)*DZ(I,K))/ &
902 (DZ(I,K)+DZ(I,K+1)))
903 SMIN = MIN(SMIN,(f(I,K)*DZ(I,K-1)+f(I,K-1)*DZ(I,K))/ &
904 (DZ(I,K)+DZ(I,K-1)), &
905 (f(I,K)*DZ(I,K+1)+f(I,K+1)*DZ(I,K))/ &
906 (DZ(I,K)+DZ(I,K+1)))
907 END IF
908
909 IF(SMIN-fn(I,K) > 0.0_SP)fn(I,K) = SMIN
910 IF(fn(I,K)-SMAX > 0.0_SP)fn(I,K) = SMAX
911
912 END DO
913 END DO nodes
914
915 WHERE(fn < 0.0_SP)fn=0.0_SP
916
917 IF(DBG_SET(DBG_SBR)) WRITE(IPT,*)"End: fct_sed"
918 End Subroutine fct_sed
919
920 !==========================================================================
921 ! Calculate Fluxes for Vertical Advection Equation
922 ! n: number of cells
923 ! c: scalar variable (1:n)
924 ! w: velocity field at cell interfaces (1:n+1)
925 ! note: we dont use face normals to construct inflow/outflow
926 ! thus we add dissipation term instead of subtracting because
927 ! positive velocity is up while computational coordinates increase
928 ! down towards bottom.
929 !==========================================================================
930 Subroutine Calc_VFlux(n,c,w,flux)
931 use mod_prec
932 implicit none
933 integer , intent(in ) :: n
934 real(sp), intent(in ) :: c(n)
935 real(sp), intent(in ) :: w(n+1)
936 real(sp), intent(out) :: flux(n+1)
937 real(sp) :: conv(n+1),diss(n+1)
938 real(sp) :: cin(-1:n+2)
939 real(sp) :: dis4
940 integer :: i
941
942 !transfer to working array
943 cin(1:n) = c(1:n)
944
945 !surface bcs (no flux)
946 cin(0) = -cin(1)
947 cin(-1) = -cin(2)
948
949 !bottom bcs (no flux)
950 cin(n+1) = -cin(n)
951 cin(n+2) = -cin(n-1)
952
953 !flux computation
954 do i=1,n+1
955 dis4 = .5*abs(w(i))
956 conv(i) = w(i)*(cin(i)+cin(i-1))/2.
957 diss(i) = dis4*(cin(i)-cin(i-1)-lim(cin(i+1)-cin(i),cin(i-1)-cin(i-2)))
958 flux(i) = conv(i)+diss(i)
959 end do
960
961 End Subroutine Calc_VFlux
962
963 !==========================================================================
964 ! Calculate LED Limiter L(u,v)
965 !==========================================================================
966 Function Lim(a,b)
967 use mod_prec
968 real(sp) lim,a,b
969 real(sp) q,R
970 real(sp) eps
971 eps = epsilon(eps)
972
973 q = 0. !1st order
974 q = 1. !minmod
975 q = 2. !van leer
976
977 R = abs( (a-b)/(abs(a)+abs(b)+eps) )**q
978 lim = .5*(1-R)*(a+b)
979
980 End Function Lim
981
982
983 End Module Scalar
为了更快速的了解计算过程，自己设置了一个只有7个节点、6个单元的简单地形，如下图所示，然后通过 printf 的方法快速了解每个变量含义。
有限体积法离散控制方程（理论）
首先介绍 FVCOM 控制方程离散，虽然其也是按照有限体积方法思想将积分降维，但是只是在平面二维方向上使用有限体积法，垂向计算使用的是有限差分法。
控制方程
控制方程离散过程
最终更新标量值Ci的程序为： fn(i,k)=(f(i,k)-xflux(i,k)/art1(i)*(deltat/dt(i)))*(dt(i)/dtfa(i))
这里deltat表示时间步长，dt为上个时间步水深，dtfa为新计算水深，这里之所以需要除以水深是因为垂向梯度项计算需要。
数值求解（结合程序分析）
首先说一下 FVCOM 选取控制体的方法，如下图所示。例如节点1所在控制体，就是由下边4条红线和上边2条黑线所组成的6面体，而节点6则是由周围12条红线组成。
FVCOM计算时候也不是按照控制体进行循环，而是按照控制边的个数循环，也就是说，像节点1和节点6这种相邻控制体，两条红色邻边各只计算一次，控制边的通量在节点1和6控制体内是大小相同、符号相反的。
1. 水平对流项计算
i1=ntrg(i) dtij(i)=dt1(i1) do k=1,kbm1 uvn(i,k) = v(i1,k)*dltxe(i) - u(i1,k)*dltye(i) end do exflux=-un*dtij(i)* ((1.0_sp+sign(1.0_sp,un))*fij2+(1.0_sp-sign(1.0_sp,un))*fij1)*0.5_sp+fxx+fyy
这里前面一项 -un*dtij(i)*((1.0_sp+sign(1.0_sp,un))*fij2+(1.0_sp-sign(1.0_sp,un))*fij1)*0.5_sp 毫无疑问就是水平对流项，((1.0_sp+sign(1.0_sp,un))*fij2+(1.0_sp-sign(1.0_sp,un))*fij1)*0.5_sp 表示其采用的是迎风格式，控制边上 fij1 及 fij2 计算采用泰勒公式达到2阶精度，泰勒公式中的一阶偏导计算在后面统一介绍
2. 水平扩散项计算
txx=0.5_sp*(pfpxd(ia)+pfpxd(ib))*viscof
tyy=0.5_sp*(pfpyd(ia)+pfpyd(ib))*viscof fxx=-dtij(i)*txx*dltye(i)
fyy= dtij(i)*tyy*dltxe(i) exflux=-un*dtij(i)* &
((1.0_sp+sign(1.0_sp,un))*fij2+(1.0_sp-sign(1.0_sp,un))*fij1)*0.5_sp+fxx+fyy
其中 viscof 为水平扩散系数，一阶偏导采用控制边相邻两个节点平均值。 fxx 和 fyy 中还包含了水深dtij，后面会在计算流量时除去。
exflux=-un*dtij(i)* &
((1.0_sp+sign(1.0_sp,un))*fij2+(1.0_sp-sign(1.0_sp,un))*fij1)*0.5_sp+fxx+fyy xflux(ia,k)=xflux(ia,k)+exflux
xflux(ib,k)=xflux(ib,k)-exflux fn(i,k)=(f(i,k)-xflux(i,k)/art1(i)*(deltat/dt(i)))*(dt(i)/dtfa(i))
3. 一阶偏导项计算
一阶偏导计算也采用格林公式将面积分降维化为线积分计算，但是平面积分所采用的控制体和上面不同，这里采用和节点相邻的所有三角形单元计算，对于边界点来说，只取相邻的两个单元。
对应的程序如下
pfpx(i) = pfpx(i) +ff1*(vy(i1)-vy(i2))
pfpy(i) = pfpy(i) +ff1*(vx(i2)-vx(i1))
pfpxd(i)= pfpxd(i)+ffd*(vy(i1)-vy(i2))
pfpyd(i)= pfpyd(i)+ffd*(vx(i2)-vx(i1)) …… PFPX(I)=PFPX(I)/ART2(I)
PFPY(I)=PFPY(I)/ART2(I)
PFPXD(I)=PFPXD(I)/ART2(I)
PFPYD(I)=PFPYD(I)/ART2(I)
对于水平对流项中控制边上Ci的二阶精度计算，只需按照泰勒公式计算即可
fij1=f(ia,k)+dxa*pfpx(ia)+dya*pfpy(ia)
fij2=f(ib,k)+dxb*pfpx(ib)+dyb*pfpy(ib)