SUNTANS 及 FVCOM 对流扩散方程求解简介［TBC］
最近接到一个任务，就是解决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)
