SUNTANS|SUNTANS 及 FVCOM 对流扩散方程求解简介［TBC］ SUNTANS及FVCOM对流扩散方程求解

最近接到一个任务，就是解决FVCOM中对流扩散计算不守衡问题。导师认为是其求解时候水平和垂向计算分开求解所导致的，目前我也没搞清到底有什么问题，反正就是让把SUNTANS的对流扩散计算挪到FVCOM中，下面就把 SUNTANS 和 FVCOM 数值求解的过程贴出来，备忘
SUNTANS模型求解过程 【SUNTANS|SUNTANS 及 FVCOM 对流扩散方程求解简介［TBC］】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, 36REAL kappa, REAL kappaH, REAL **kappa_tv, REAL theta, 37REAL **src1, REAL **src2, REAL *Ftop, REAL *Fbot, int alpha_top, int alpha_bot, 38MPI_Comm comm, int myproc, int checkflag, int TVDscheme) 39 { 40int i, iptr, j, jptr, ib, k, nf, ktop; 41int Nc=grid->Nc, normal, nc1, nc2, ne; 42REAL df, dg, Ac, dt=prop->dt, fab, *a, *b, *c, *d, *ap, *am, *bd, *uflux, dznew, mass, *sp, *temp; 43REAL smin, smax, div_local, div_da; 44int k1, k2, kmin, imin, kmax, imax, mincount, maxcount, allmincount, allmaxcount, flag; 45 46prop->TVD = TVDscheme; 47// These are used mostly debugging to turn on/off vertical and horizontal TVD. 48prop->horiTVD = 1; 49prop->vertTVD = 1; 50 51ap = phys->ap; 52am = phys->am; 53bd = phys->bp; 54temp = phys->bm; 55a = phys->a; 56b = phys->b; 57c = phys->c; 58d = phys->d; 59 60// Never use AB2 61if(1) { 62fab=1; 63for(i=0; iNc; i++) 64for(k=0; kNk[i]; k++) 65Cn[i][k]=0; 66} else 67fab=1.5; 68 69for(i=0; iNk[i]; k++) 71phys->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]; iptrcelldist[1]; iptr++) { 76for(iptr=grid->celldist[0]; iptrcelldist[2]; iptr++) { 77i = grid->cellp[iptr]; 78 79for(k=grid->ctop[i]; kNk[i]; k++) 80phys->stmp2[i][k]=0; 81} 82 83if(boundary_scal) { 84for(jptr=grid->edgedist[2]; jptredgedist[5]; jptr++) { 85j = grid->edgep[jptr]; 86 87ib = 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. 91for(k=grid->ctop[ib]; kNk[ib]; k++) 92phys->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 98if(prop->TVD && prop->horiTVD) 99HorizontalFaceScalars(grid,phys,prop,scal,boundary_scal,prop->TVD,comm,myproc); 100 101//for(iptr=grid->celldist[0]; iptrcelldist[1]; iptr++) { 102for(iptr=grid->celldist[0]; iptrcelldist[2]; iptr++) { 103i = grid->cellp[iptr]; 104Ac = grid->Ac[i]; 105 106if(grid->ctop[i]>=grid->ctopold[i]) { 107ktop=grid->ctop[i]; 108dznew=grid->dzz[i][ktop]; 109} else { 110ktop=grid->ctopold[i]; 111dznew=0; 112for(k=grid->ctop[i]; k<=grid->ctopold[i]; k++) 113dznew+=grid->dzz[i][k]; 114} 115 116// These are the advective components of the tridiagonal 117// at the new time step. 118if(!(prop->TVD && prop->vertTVD)) 119for(k=0; kNk[i]+1; k++) { 120ap[k] = 0.5*(wnew[i][k]+fabs(wnew[i][k])); 121am[k] = 0.5*(wnew[i][k]-fabs(wnew[i][k])); 122} 123else// Compute the ap/am for TVD schemes 124GetApAm(ap,am,phys->wp,phys->wm,phys->Cp,phys->Cm,phys->rp,phys->rm, 125wnew,grid->dzz,scal,i,grid->Nk[i],ktop,prop->dt,prop->TVD); 126 127for(k=ktop+1; kNk[i]; k++) { 128a[k-ktop]=theta*dt*am[k]; 129b[k-ktop]=grid->dzz[i][k]+theta*dt*(ap[k]-am[k+1]); 130c[k-ktop]=-theta*dt*ap[k+1]; 131} 132 133// Top cell advection 134a[0]=0; 135b[0]=dznew-theta*dt*am[ktop+1]; 136c[0]=-theta*dt*ap[ktop+1]; 137 138// Bottom cell no-flux boundary condition for advection 139b[(grid->Nk[i]-1)-ktop]+=c[(grid->Nk[i]-1)-ktop]; 140 141// Implicit vertical diffusion terms 142for(k=ktop+1; kNk[i]; k++) 143bd[k]=(2.0*kappa+kappa_tv[i][k-1]+kappa_tv[i][k])/ 144(grid->dzz[i][k-1]+grid->dzz[i][k]); 145 146for(k=ktop+1; kNk[i]-1; k++) { 147a[k-ktop]-=theta*dt*bd[k]; 148b[k-ktop]+=theta*dt*(bd[k]+bd[k+1]); 149c[k-ktop]-=theta*dt*bd[k+1]; 150} 151if(src1) 152for(k=ktop; kNk[i]; k++) 153b[k-ktop]+=grid->dzz[i][k]*src1[i][k]*theta*dt; 154 155// Diffusive fluxes only when more than 1 layer 156if(ktopNk[i]-1) { 157// Top cell diffusion 158b[0]+=theta*dt*(bd[ktop+1]+2*alpha_top*bd[ktop+1]); 159c[0]-=theta*dt*bd[ktop+1]; 160 161// Bottom cell diffusion 162a[(grid->Nk[i]-1)-ktop]-=theta*dt*bd[grid->Nk[i]-1]; 163b[(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[] 167for(k=ktop+1; kNk[i]; k++) 168d[k-ktop]=grid->dzzold[i][k]*phys->stmp[i][k]; 169if(src1) 170for(k=ktop+1; kNk[i]; k++) 171d[k-ktop]-=src1[i][k]*(1-theta)*dt*grid->dzzold[i][k]*phys->stmp[i][k]; 172 173d[0]=0; 174if(grid->ctopold[i]<=grid->ctop[i]) { 175for(k=grid->ctopold[i]; k<=grid->ctop[i]; k++) 176d[0]+=grid->dzzold[i][k]*phys->stmp[i][k]; 177if(src1) 178for(k=grid->ctopold[i]; k<=grid->ctop[i]; k++) 179d[0]-=src1[i][k]*(1-theta)*dt*grid->dzzold[i][k]*phys->stmp[i][k]; 180} else { 181d[0]=grid->dzzold[i][ktop]*phys->stmp[i][ktop]; 182if(src1) 183d[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 188if(!(prop->TVD && prop->vertTVD)) 189for(k=0; kNk[i]+1; k++) { 190ap[k] = 0.5*(phys->wtmp2[i][k]+fabs(phys->wtmp2[i][k])); 191am[k] = 0.5*(phys->wtmp2[i][k]-fabs(phys->wtmp2[i][k])); 192} 193else // Compute the ap/am for TVD schemes 194GetApAm(ap,am,phys->wp,phys->wm,phys->Cp,phys->Cm,phys->rp,phys->rm, 195phys->wtmp2,grid->dzzold,phys->stmp,i,grid->Nk[i],ktop,prop->dt,prop->TVD); 196 197// Explicit advection and diffusion 198for(k=ktop+1; kNk[i]-1; k++) 199d[k-ktop]-=(1-theta)*dt*(am[k]*phys->stmp[i][k-1]+ 200(ap[k]-am[k+1])*phys->stmp[i][k]- 201ap[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 206if(ktopNk[i]-1) { 207//Flux through bottom of top cell 208k=ktop; 209d[0]=d[0]-(1-theta)*dt*(-am[k+1]*phys->stmp[i][k]- 210ap[k+1]*phys->stmp[i][k+1])+ 211(1-theta)*dt*(-(2*alpha_top*bd[k+1]+bd[k+1])*phys->stmp[i][k]+ 212bd[k+1]*phys->stmp[i][k+1]); 213if(Ftop) d[0]+=dt*(1-alpha_top+2*alpha_top*bd[k+1])*Ftop[i]; 214 215// Through top of bottom cell 216k=grid->Nk[i]-1; 217d[k-ktop]-=(1-theta)*dt*(am[k]*phys->stmp[i][k-1]+ 218ap[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]); 221if(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. 225if(grid->ctop[i]<=grid->ctopold[i]) { 226for(k=grid->ctop[i]; k<=grid->ctopold[i]; k++) 227d[0]+=(1-fab)*Cn[i][grid->ctopold[i]]/(1+fabs(grid->ctop[i]-grid->ctopold[i])); 228for(k=grid->ctopold[i]+1; kNk[i]; k++) 229d[k-grid->ctopold[i]]+=(1-fab)*Cn[i][k]; 230} else { 231for(k=grid->ctopold[i]; k<=grid->ctop[i]; k++) 232d[0]+=(1-fab)*Cn[i][k]; 233for(k=grid->ctop[i]+1; kNk[i]; k++) 234d[k-grid->ctop[i]]+=(1-fab)*Cn[i][k]; 235} 236 237for(k=0; kctop[i]; k++) 238Cn[i][k]=0; 239 240if(src2) 241for(k=grid->ctop[i]; kNk[i]; k++) 242Cn[i][k-ktop]=dt*src2[i][k]*grid->dzzold[i][k]; 243else 244for(k=grid->ctop[i]; kNk[i]; k++) 245Cn[i][k]=0; 246 247// Now create the source term for the current time step 248for(k=0; kNk[i]; k++) 249ap[k]=0; 250 251for(nf=0; nfnfaces[i]; nf++) { 252ne = grid->face[i*grid->maxfaces+nf]; 253normal = grid->normal[i*grid->maxfaces+nf]; 254df = grid->df[ne]; 255dg = grid->dg[ne]; 256nc1 = grid->grad[2*ne]; 257nc2 = grid->grad[2*ne+1]; 258if(nc1==-1) nc1=nc2; 259if(nc2==-1) { 260nc2=nc1; 261if(boundary_scal && (grid->mark[ne]==2 || grid->mark[ne]==3)) 262sp=phys->stmp2[nc1]; 263else 264sp=phys->stmp[nc1]; 265} else 266sp=phys->stmp[nc2]; 267 268if(!(prop->TVD && prop->horiTVD)) { 269for(k=0; kNke[ne]; k++) 270temp[k]=UpWind(phys->utmp2[ne][k], 271phys->stmp[nc1][k], 272sp[k]); 273} else { 274for(k=0; kNke[ne]; k++) 275if(phys->utmp2[ne][k]>0) 276temp[k]=phys->SfHp[ne][k]; 277else 278temp[k]=phys->SfHm[ne][k]; 279} 280 281for(k=0; kNke[ne]; k++) 282ap[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 286for(k=ktop+1; kNk[i]; k++) 287Cn[i][k-ktop]-=ap[k]; 288 289for(k=0; k<=ktop; k++) 290Cn[i][0]-=ap[k]; 291 292// Add on the source from the current time step to the rhs. 293for(k=0; kNk[i]-ktop; k++) 294d[k]+=fab*Cn[i][k]; 295 296// Add on the volume correction if h was < -d 297/* 298if(grid->ctop[i]==grid->Nk[i]-1) 299d[grid->Nk[i]-ktop-1]+=phys->hcorr[i]*phys->stmp[i][grid->ctop[i]]; 300*/ 301 302for(k=ktop; kNk[i]; k++) 303ap[k]=Cn[i][k-ktop]; 304for(k=0; k<=ktop; k++) 305Cn[i][k]=0; 306for(k=ktop+1; kNk[i]; k++) 307Cn[i][k]=ap[k]; 308for(k=grid->ctop[i]; k<=ktop; k++) 309Cn[i][k]=ap[ktop]/(1+fabs(grid->ctop[i]-ktop)); 310 311if(grid->Nk[i]-ktop>1) 312TriSolve(a,b,c,d,&(scal[i][ktop]),grid->Nk[i]-ktop); 313else if(prop->n>1) { 314if(b[0]>0 && phys->active[i]) 315scal[i][ktop]=d[0]/b[0]; 316else 317scal[i][ktop]=0; 318} 319 320for(k=0; kctop[i]; k++) 321scal[i][k]=0; 322 323for(k=grid->ctop[i]; kctopold[i]; k++) 324scal[i][k]=scal[i][ktop]; 325} 326 327// Code to check divergence change CHECKCONSISTENCY to 1 in suntans.h 328if(CHECKCONSISTENCY && checkflag) { 329 330if(prop->n==1+prop->nstart) { 331smin=INFTY; 332smax=-INFTY; 333for(i=0; iNc; i++) { 334for(k=grid->ctop[i]; kNk[i]; k++) { 335if(phys->stmp[i][k]>smax) { 336smax=phys->stmp[i][k]; 337imax=i; 338kmax=k; 339} 340if(phys->stmp[i][k]stmp[i][k]; 342imin=i; 343kmin=k; 344} 345} 346} 347MPI_Reduce(&smin,&smin_value,1,MPI_DOUBLE,MPI_MIN,0,comm); 348MPI_Reduce(&smax,&smax_value,1,MPI_DOUBLE,MPI_MAX,0,comm); 349MPI_Bcast(&smin_value,1,MPI_DOUBLE,0,comm); 350MPI_Bcast(&smax_value,1,MPI_DOUBLE,0,comm); 351 352if(myproc==0) 353printf("Minimum scalar: %.2f, maximum: %.2f\n",smin_value,smax_value); 354} 355 356//for(iptr=grid->celldist[0]; iptrcelldist[1]; iptr++) { 357for(iptr=grid->celldist[0]; iptrcelldist[2]; iptr++) { 358i = grid->cellp[iptr]; 359 360flag=0; 361for(nf=0; nfnfaces[i]; nf++) { 362if(grid->mark[grid->face[i*grid->maxfaces+nf]]==2 || 363grid->mark[grid->face[i*grid->maxfaces+nf]]==3) { 364flag=1; 365break; 366} 367} 368 369if(!flag) { 370div_da=0; 371 372for(k=0; kNk[i]; k++) { 373div_da+=grid->Ac[i]*(grid->dzz[i][k]-grid->dzzold[i][k])/prop->dt; 374 375div_local=0; 376for(nf=0; nfnfaces[i]; nf++) { 377ne=grid->face[i*grid->maxfaces+nf]; 378div_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} 381div_da+=div_local; 382div_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 385if(k>=grid->ctop[i]) { 386if(fabs(div_local)>SMALL_CONSISTENCY && grid->dzz[imin][0]>DRYCELLHEIGHT) 387printf("Step: %d, proc: %d, locally-divergent at %d, %d, div=%e\n", 388prop->n,myproc,i,k,div_local); 389} 390} 391if(fabs(div_da)>SMALL_CONSISTENCY && phys->h[i]+grid->dv[i]>DRYCELLHEIGHT) 392printf("Step: %d, proc: %d, Depth-Ave divergent at i=%d, div=%e\n", 393prop->n,myproc,i,div_da); 394} 395} 396 397mincount=0; 398maxcount=0; 399smin=INFTY; 400smax=-INFTY; 401//for(iptr=grid->celldist[0]; iptrcelldist[1]; iptr++) { 402for(iptr=grid->celldist[0]; iptrcelldist[2]; iptr++) { 403i = grid->cellp[iptr]; 404 405flag=0; 406for(nf=0; nfnfaces[i]; nf++) { 407if(grid->mark[grid->face[i*grid->maxfaces+nf]]==2 || grid->mark[grid->face[i*grid->maxfaces+nf]]==3) { 408flag=1; 409break; 410} 411} 412 413if(!flag) { 414for(k=grid->ctop[i]; kNk[i]; k++) { 415if(scal[i][k]>smax) { 416smax=scal[i][k]; 417imax=i; 418kmax=k; 419} 420if(scal[i][k]smax_value+SMALL_CONSISTENCY && grid->dzz[i][k]>DRYCELLHEIGHT) 427maxcount++; 428if(scal[i][k]dzz[i][k]>DRYCELLHEIGHT) 429mincount++; 430} 431} 432} 433MPI_Reduce(&mincount,&allmincount,1,MPI_INT,MPI_SUM,0,comm); 434MPI_Reduce(&maxcount,&allmaxcount,1,MPI_INT,MPI_SUM,0,comm); 435 436if(mincount!=0 || maxcount!=0) 437printf("Not CWC, step: %d, proc: %d, smin = %e at i=%d,H=%e, smax = %e at i=%d,H=%e\n", 438prop->n,myproc, 439smin,imin,phys->h[imin]+grid->dv[imin], 440smax,imax,phys->h[imax]+grid->dv[imax]); 441 442if(myproc==0 && (allmincount !=0 || allmaxcount !=0)) 443printf("Total number of CWC violations (all procs): ss_max: %d\n", 444allmincount,allmaxcount); 445} 446}

View Code

下面介绍解读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 ORIMPLIED WARRANTIES, INCLUDING,BUT NOTLIMITED TO, 19 ! THE IMPLIED WARRANTIES OF MERCHANTABILITY ANDFITNESS 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 40logical, parameter :: debug = .true. 41 42contains 43 !==============================================================================| 44 ! Calculate Horizontal Advection and Diffusion For Scalar (f)| 45 !==============================================================================| 46Subroutine Adv_Scal(f,fn,d_fdis,fdis,d_fflux,fflux_obc,deltat,source) 47 !------------------------------------------------------------------------------| 48 49use all_vars 50use lims, only: m,mt,n,nt,kbm1,kb 51use bcs 52use mod_obcs 53 # if defined (MULTIPROCESSOR) 54use mod_par 55 # endif 56 # if defined (WET_DRY) 57use mod_wd 58 # endif 59 60 # if defined (THIN_DAM) 61use mod_dam, only : kdam,N_DAM_MATCH,IS_DAM 62 # endif 63 64implicit none 65real(sp), intent(in ), dimension(0:mt,kb):: f 66real(sp), intent(out), dimension(0:mt,kb):: fn 67integer , intent(in ):: d_fdis 68real(sp), intent(in ), dimension(d_fdis):: fdis 69integer , intent(in ):: d_fflux 70real(sp), intent(out), dimension(d_fflux,kbm1) :: fflux_obc 71real(sp), intent(in ):: deltat 72logical , intent(in ):: source 73 74!----------------local-------------------------------------- 75real(sp), dimension(0:mt,kb):: xflux,xflux_adv 76real(sp), dimension(m):: pupx,pupy,pvpx,pvpy 77real(sp), dimension(m):: pfpx,pfpy,pfpxd,pfpyd,viscoff 78real(sp), dimension(3*nt):: dtij 79real(sp), dimension(3*nt,kbm1) :: uvn 80real(sp), dimension(kb):: vflux 81real(sp) :: utmp,vtmp,sitai,ffd,ff1,x11,y11,x22,y22,x33,y33 82real(sp) :: tmp1,tmp2,xi,yi 83real(sp) :: dxa,dya,dxb,dyb,fij1,fij2,un 84real(sp) :: txx,tyy,fxx,fyy,viscof,exflux,temp,fpoint 85real(sp) :: fact,fm1,fmean 86integer:: i,i1,i2,ia,ib,j,j1,j2,k,jtmp,jj 87 # if defined (SPHERICAL) 88real(sp) :: ty,txpi,typi 89 # endif 90 91 # if defined (THIN_DAM) 92INTEGER:: NX 93real(sp) :: tmpflx 94real(sp),dimension(kb) :: wvel 95 # endif 96 97 98 !------------------------------------------------------------------------------! 99 100 !------------------------------------------------------- 101 !Calculate Mean Values 102 !------------------------------------------------------- 103 104fmean = sum(f(1:m,1:kbm1))/float(m*kbm1) 105 106 !------------------------------------------------------- 107 !Initialize Multipliers to Control Horizontal Diff 108 !------------------------------------------------------- 109 110fact = 0.0_sp 111fm1= 1.0_sp 112if(HORIZONTAL_MIXING_TYPE == 'closure') then 113fact = 1.0_sp 114fm1= 0.0_sp 115end if 116 117 !------------------------------------------------------- 118 !Initialize Fluxes 119 !------------------------------------------------------- 120xflux= 0.0_sp 121xflux_adv = 0.0_sp 122 123 !------------------------------------------------------- 124 !Calculate Normal Velocity on Control Volume Edges 125 !------------------------------------------------------- 126 !!# if !defined (WET_DRY) 127do i=1,ncv 128i1=ntrg(i) 129dtij(i)=dt1(i1) 130do k=1,kbm1 131uvn(i,k) = v(i1,k)*dltxe(i) - u(i1,k)*dltye(i) 132 133 #if defined(PLBC) 134uvn(i,k) =- u(i1,k)*dltye(i) 135 #endif 136 137end do 138end 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 153do k=1,kbm1 154pfpx= 0.0_sp 155pfpy= 0.0_sp 156pfpxd = 0.0_sp 157pfpyd = 0.0_sp 158do i=1,m 159do j=1,ntsn(i)-1 160i1=nbsn(i,j) 161i2=nbsn(i,j+1) 162 163 #if defined (WET_DRY) 164IF(ISWETN(I1) == 0 .AND. ISWETN(I2) == 1)THEN 165FFD=0.5_SP*(f(I,K)+f(I2,K)) 166FF1=0.5_SP*(f(I,K)+f(I2,K)) 167ELSE IF(ISWETN(I1) == 1 .AND. ISWETN(I2) == 0)THEN 168FFD=0.5_SP*(f(I1,K)+f(I,K)) 169FF1=0.5_SP*(f(I1,K)+f(I,K)) 170ELSE IF(ISWETN(I1) == 0 .AND. ISWETN(I2) == 0)THEN 171FFD=0.5_SP*(f(I,K)+f(I,K)) 172FF1=0.5_SP*(f(I,K)+f(I,K)) 173ELSE 174FFD=0.5_SP*(f(I1,K)+f(I2,K)) 175FF1=0.5_SP*(f(I1,K)+f(I2,K)) 176END IF 177 #else 178ffd=0.5_sp*(f(i1,k)+f(i2,k)) !-fmean1(i1,k)-fmean1(i2,k)) 179ff1=0.5_sp*(f(i1,k)+f(i2,k)) 180 #endif 181 182 #if defined (SPHERICAL) 183ty=0.5_sp*(vy(i1)+vy(i2)) 184txpi=(vx(i2)-vx(i1))*tpi*cos(deg2rad*ty) 185typi=(vy(i1)-vy(i2))*tpi 186pfpx(i)=pfpx(i)+ff1*typi 187pfpy(i)=pfpy(i)+ff1*txpi 188pfpxd(i)=pfpxd(i)+ffd*typi 189pfpyd(i)=pfpyd(i)+ffd*txpi 190 #else 191pfpx(i) = pfpx(i) +ff1*(vy(i1)-vy(i2)) 192pfpy(i) = pfpy(i) +ff1*(vx(i2)-vx(i1)) 193pfpxd(i)= pfpxd(i)+ffd*(vy(i1)-vy(i2)) 194pfpyd(i)= pfpyd(i)+ffd*(vx(i2)-vx(i1)) 195 #endif 196end do 197 198 ! gather all neighboring control volumes connecting at dam node 199 # if defined (THIN_DAM) 200IF(IS_DAM(I)==1.AND.K<=KDAM(I))THEN 201DO NX=1,N_DAM_MATCH(I,1) 202DO J=1,NTSN(N_DAM_MATCH(I,NX+1))-1 203I1=NBSN(N_DAM_MATCH(I,NX+1),J) 204I2=NBSN(N_DAM_MATCH(I,NX+1),J+1) 205FFD=0.5_SP*(f(I1,K)+f(I2,K)) !-SMEAN1(I1,K)-SMEAN1(I2,K)) 206FF1=0.5_SP*(f(I1,K)+f(I2,K)) 207 #if defined (SPHERICAL) 208ty=0.5_sp*(vy(i1)+vy(i2)) 209txpi=(vx(i2)-vx(i1))*tpi*cos(deg2rad*ty) 210typi=(vy(i1)-vy(i2))*tpi 211pfpx(i)=pfpx(i)+ff1*typi 212pfpy(i)=pfpy(i)+ff1*txpi 213pfpxd(i)=pfpxd(i)+ffd*typi 214pfpyd(i)=pfpyd(i)+ffd*txpi 215 #else 216pfpx(i) = pfpx(i) +ff1*(vy(i1)-vy(i2)) 217pfpy(i) = pfpy(i) +ff1*(vx(i2)-vx(i1)) 218pfpxd(i)= pfpxd(i)+ffd*(vy(i1)-vy(i2)) 219pfpyd(i)= pfpyd(i)+ffd*(vx(i2)-vx(i1)) 220 #endif 221END DO 222END DO 223END IF 224 # endif 225 226 # if !defined (THIN_DAM) 227pfpx(i)=pfpx(i )/art2(i) 228pfpy(i)=pfpy(i )/art2(i) 229pfpxd(i) =pfpxd(i)/art2(i) 230pfpyd(i) =pfpyd(i)/art2(i) 231 # else 232IF(IS_DAM(I)==1.AND.K<=KDAM(I))THEN 233PFPX(I)=PFPX(I)/(ART2(I)+SUM(ART2(N_DAM_MATCH(I,2:1+N_DAM_MATCH(I,1))))) 234PFPY(I)=PFPY(I)/(ART2(I)+SUM(ART2(N_DAM_MATCH(I,2:1+N_DAM_MATCH(I,1))))) 235PFPXD(I)=PFPXD(I)/(ART2(I)+SUM(ART2(N_DAM_MATCH(I,2:1+N_DAM_MATCH(I,1))))) 236PFPYD(I)=PFPYD(I)/(ART2(I)+SUM(ART2(N_DAM_MATCH(I,2:1+N_DAM_MATCH(I,1))))) 237ELSE 238PFPX(I)=PFPX(I)/ART2(I) 239PFPY(I)=PFPY(I)/ART2(I) 240PFPXD(I)=PFPXD(I)/ART2(I) 241PFPYD(I)=PFPYD(I)/ART2(I) 242END IF 243 # endif 244 245end do 246 247if(k == kbm1)then 248do i=1,m 249pfpxb(i) = pfpx(i) 250pfpyb(i) = pfpy(i) 251end do 252end if 253 254do i=1,m 255pupx(i)=0.0_sp 256pupy(i)=0.0_sp 257pvpx(i)=0.0_sp 258pvpy(i)=0.0_sp 259j=1 260i1=nbve(i,j) 261jtmp=nbvt(i,j) 262j1=jtmp+1-(jtmp+1)/4*3 263j2=jtmp+2-(jtmp+2)/4*3 264x11=0.5_sp*(vx(i)+vx(nv(i1,j1))) 265y11=0.5_sp*(vy(i)+vy(nv(i1,j1))) 266x22=xc(i1) 267y22=yc(i1) 268x33=0.5_sp*(vx(i)+vx(nv(i1,j2))) 269y33=0.5_sp*(vy(i)+vy(nv(i1,j2))) 270 271 #if defined (SPHERICAL) 272ty=0.5_sp*(y11+y33) 273txpi=(x33-x11)*tpi*cos(deg2rad*ty) 274typi=(y11-y33)*tpi 275pupx(i)=pupx(i)+u(i1,k)*typi 276pupy(i)=pupy(i)+u(i1,k)*txpi 277pvpx(i)=pvpx(i)+v(i1,k)*typi 278pvpy(i)=pvpy(i)+v(i1,k)*txpi 279 #else 280pupx(i)=pupx(i)+u(i1,k)*(y11-y33) 281pupy(i)=pupy(i)+u(i1,k)*(x33-x11) 282pvpx(i)=pvpx(i)+v(i1,k)*(y11-y33) 283pvpy(i)=pvpy(i)+v(i1,k)*(x33-x11) 284 #endif 285 286if(isonb(i) /= 0) then 287 #if defined (SPHERICAL) 288ty=0.5_sp*(vy(i)+y11) 289txpi=(x11-vx(i))*tpi*cos(deg2rad*ty) 290typi=(vy(i)-y11)*tpi 291pupx(i)=pupx(i)+u(i1,k)*typi 292pupy(i)=pupy(i)+u(i1,k)*txpi 293pvpx(i)=pvpx(i)+v(i1,k)*typi 294pvpy(i)=pvpy(i)+v(i1,k)*txpi 295 #else 296pupx(i)=pupx(i)+u(i1,k)*(vy(i)-y11) 297pupy(i)=pupy(i)+u(i1,k)*(x11-vx(i)) 298pvpx(i)=pvpx(i)+v(i1,k)*(vy(i)-y11) 299pvpy(i)=pvpy(i)+v(i1,k)*(x11-vx(i)) 300 #endif 301end if 302 303do j=2,ntve(i)-1 304i1=nbve(i,j) 305jtmp=nbvt(i,j) 306j1=jtmp+1-(jtmp+1)/4*3 307j2=jtmp+2-(jtmp+2)/4*3 308x11=0.5_sp*(vx(i)+vx(nv(i1,j1))) 309y11=0.5_sp*(vy(i)+vy(nv(i1,j1))) 310x22=xc(i1) 311y22=yc(i1) 312x33=0.5_sp*(vx(i)+vx(nv(i1,j2))) 313y33=0.5_sp*(vy(i)+vy(nv(i1,j2))) 314 315 #if defined (SPHERICAL) 316ty=0.5_sp*(y11+y33) 317txpi=(x33-x11)*tpi*COS(deg2rad*TY) 318typi=(y11-y33)*tpi 319pupx(i)=pupx(i)+u(i1,k)*typi 320pupy(i)=pupy(i)+u(i1,k)*txpi 321pvpx(i)=pvpx(i)+v(i1,k)*typi 322pvpy(i)=pvpy(i)+v(i1,k)*txpi 323 #else 324pupx(i)=pupx(i)+u(i1,k)*(y11-y33) 325pupy(i)=pupy(i)+u(i1,k)*(x33-x11) 326pvpx(i)=pvpx(i)+v(i1,k)*(y11-y33) 327pvpy(i)=pvpy(i)+v(i1,k)*(x33-x11) 328 #endif 329end do 330j=ntve(i) 331i1=nbve(i,j) 332jtmp=nbvt(i,j) 333j1=jtmp+1-(jtmp+1)/4*3 334j2=jtmp+2-(jtmp+2)/4*3 335x11=0.5_sp*(vx(i)+vx(nv(i1,j1))) 336y11=0.5_sp*(vy(i)+vy(nv(i1,j1))) 337x22=xc(i1) 338y22=yc(i1) 339x33=0.5_sp*(vx(i)+vx(nv(i1,j2))) 340y33=0.5_sp*(vy(i)+vy(nv(i1,j2))) 341 342 #if defined (SPHERICAL) 343ty=0.5*(Y11+Y33) 344txpi=(x33-x11)*tpi*cos(deg2rad*TY) 345typi=(y11-y33)*tpi 346pupx(i)=pupx(i)+u(i1,k)*typi 347pupy(i)=pupy(i)+u(i1,k)*txpi 348pvpx(i)=pvpx(i)+v(i1,k)*typi 349pvpy(i)=pvpy(i)+v(i1,k)*txpi 350 #else 351pupx(i)=pupx(i)+u(i1,k)*(y11-y33) 352pupy(i)=pupy(i)+u(i1,k)*(x33-x11) 353pvpx(i)=pvpx(i)+v(i1,k)*(y11-y33) 354pvpy(i)=pvpy(i)+v(i1,k)*(x33-x11) 355 #endif 356 357if(isonb(i) /= 0) then 358 #if defined (SPHERICAL) 359ty=0.5*(Y11+VY(I)) 360txpi=(VX(I)-X11)*tpi*COS(deg2rad*ty) 361typi=(Y11-VY(I))*tpi 362pupx(i)=pupx(i)+u(i1,k)*typi 363pupy(i)=pupy(i)+u(i1,k)*txpi 364pvpx(i)=pvpx(i)+v(i1,k)*typi 365pvpy(i)=pvpy(i)+v(i1,k)*txpi 366 #else 367pupx(i)=pupx(i)+u(i1,k)*(y11-vy(i)) 368pupy(i)=pupy(i)+u(i1,k)*(vx(i)-x11) 369pvpx(i)=pvpx(i)+v(i1,k)*(y11-vy(i)) 370pvpy(i)=pvpy(i)+v(i1,k)*(vx(i)-x11) 371 #endif 372end if 373pupx(i)=pupx(i)/art1(i) 374pupy(i)=pupy(i)/art1(i) 375pvpx(i)=pvpx(i)/art1(i) 376pvpy(i)=pvpy(i)/art1(i) 377tmp1=pupx(i)**2+pvpy(i)**2 378tmp2=0.5_sp*(pupy(i)+pvpx(i))**2 379viscoff(i)=sqrt(tmp1+tmp2)*art1(i) 380end do 381 !if(k == kbm1) then 382 !ah_bottom(1:m) = horcon*(fact*viscoff(1:m) + fm1) 383 !end if 384 385 386do i=1,ncv_i 387ia=niec(i,1) 388ib=niec(i,2) 389xi=0.5_sp*(xije(i,1)+xije(i,2)) 390yi=0.5_sp*(yije(i,1)+yije(i,2)) 391 #if defined (SPHERICAL) 392ty=0.5_sp*(yi+vy(ia)) 393dxa=(xi-vx(ia))*tpi*cos(deg2rad*ty) 394dya=(yi-vy(ia))*tpi 395ty=0.5*(YI+VY(IB)) 396DXB=(XI-VX(IB))*tpi*COS(deg2rad*ty) 397DYB=(YI-VY(IB))*tpi 398 #else 399dxa=xi-vx(ia) 400dya=yi-vy(ia) 401dxb=xi-vx(ib) 402dyb=yi-vy(ib) 403 #endif 404fij1=f(ia,k)+dxa*pfpx(ia)+dya*pfpy(ia) 405fij2=f(ib,k)+dxb*pfpx(ib)+dyb*pfpy(ib) 406un=uvn(i,k) 407 408 !viscof=horcon*(fact*(viscoff(ia)+viscoff(ib))*0.5_sp + fm1) 409VISCOF=(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 411txx=0.5_sp*(pfpxd(ia)+pfpxd(ib))*viscof 412tyy=0.5_sp*(pfpyd(ia)+pfpyd(ib))*viscof 413 414fxx=-dtij(i)*txx*dltye(i) 415fyy= dtij(i)*tyy*dltxe(i) 416 417 # if defined (PLBC) 418fyy=0.0_SP 419 # endif 420 421exflux=-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 424xflux(ia,k)=xflux(ia,k)+exflux 425xflux(ib,k)=xflux(ib,k)-exflux 426 427xflux_adv(ia,k)=xflux_adv(ia,k)+(exflux-fxx-fyy) 428xflux_adv(ib,k)=xflux_adv(ib,k)-(exflux-fxx-fyy) 429 430 #if defined (THIN_DAM) 431IF(K<=KDAM(IA).AND.IS_DAM(IA)==1)THEN 432IF(N_DAM_MATCH(IA,1)==1)THEN 433XFLUX(N_DAM_MATCH(IA,2),K) = XFLUX(N_DAM_MATCH(IA,2),K) + EXFLUX 434XFLUX_ADV(N_DAM_MATCH(IA,2),K) = XFLUX_ADV(N_DAM_MATCH(IA,2),K) +(EXFLUX-FXX-FYY) 435END IF 436IF(N_DAM_MATCH(IA,1)==2)THEN 437XFLUX(N_DAM_MATCH(IA,2),K) = XFLUX(N_DAM_MATCH(IA,2),K) + EXFLUX 438XFLUX(N_DAM_MATCH(IA,3),K) = XFLUX(N_DAM_MATCH(IA,3),K) + EXFLUX 439XFLUX_ADV(N_DAM_MATCH(IA,2),K) = XFLUX_ADV(N_DAM_MATCH(IA,2),K) +(EXFLUX-FXX-FYY) 440XFLUX_ADV(N_DAM_MATCH(IA,3),K) = XFLUX_ADV(N_DAM_MATCH(IA,3),K) +(EXFLUX-FXX-FYY) 441END IF 442IF(N_DAM_MATCH(IA,1)==3)THEN 443XFLUX(N_DAM_MATCH(IA,2),K) = XFLUX(N_DAM_MATCH(IA,2),K) + EXFLUX 444XFLUX(N_DAM_MATCH(IA,3),K) = XFLUX(N_DAM_MATCH(IA,3),K) + EXFLUX 445XFLUX(N_DAM_MATCH(IA,4),K) = XFLUX(N_DAM_MATCH(IA,4),K) + EXFLUX 446XFLUX_ADV(N_DAM_MATCH(IA,2),K) = XFLUX_ADV(N_DAM_MATCH(IA,2),K) +(EXFLUX-FXX-FYY) 447XFLUX_ADV(N_DAM_MATCH(IA,3),K) = XFLUX_ADV(N_DAM_MATCH(IA,3),K) +(EXFLUX-FXX-FYY) 448XFLUX_ADV(N_DAM_MATCH(IA,4),K) = XFLUX_ADV(N_DAM_MATCH(IA,4),K) +(EXFLUX-FXX-FYY) 449END IF 450END IF 451IF(K<=KDAM(IB).AND.IS_DAM(IB)==1)THEN 452IF(N_DAM_MATCH(IB,1)==1)THEN 453XFLUX(N_DAM_MATCH(IB,2),K) = XFLUX(N_DAM_MATCH(IB,2),K) - EXFLUX 454XFLUX_ADV(N_DAM_MATCH(IB,2),K) = XFLUX_ADV(N_DAM_MATCH(IB,2),K) - (EXFLUX-FXX-FYY) 455END IF 456IF(N_DAM_MATCH(IB,1)==2)THEN 457XFLUX(N_DAM_MATCH(IB,2),K) = XFLUX(N_DAM_MATCH(IB,2),K) - EXFLUX 458XFLUX(N_DAM_MATCH(IB,3),K) = XFLUX(N_DAM_MATCH(IB,3),K) - EXFLUX 459XFLUX_ADV(N_DAM_MATCH(IB,2),K) = XFLUX_ADV(N_DAM_MATCH(IB,2),K) - (EXFLUX-FXX-FYY) 460XFLUX_ADV(N_DAM_MATCH(IB,3),K) = XFLUX_ADV(N_DAM_MATCH(IB,3),K) - (EXFLUX-FXX-FYY) 461END IF 462IF(N_DAM_MATCH(IB,1)==3)THEN 463XFLUX(N_DAM_MATCH(IB,2),K) = XFLUX(N_DAM_MATCH(IB,2),K) - EXFLUX 464XFLUX(N_DAM_MATCH(IB,3),K) = XFLUX(N_DAM_MATCH(IB,3),K) - EXFLUX 465XFLUX(N_DAM_MATCH(IB,4),K) = XFLUX(N_DAM_MATCH(IB,4),K) - EXFLUX 466XFLUX_ADV(N_DAM_MATCH(IB,2),K) = XFLUX_ADV(N_DAM_MATCH(IB,2),K) - (EXFLUX-FXX-FYY) 467XFLUX_ADV(N_DAM_MATCH(IB,3),K) = XFLUX_ADV(N_DAM_MATCH(IB,3),K) - (EXFLUX-FXX-FYY) 468XFLUX_ADV(N_DAM_MATCH(IB,4),K) = XFLUX_ADV(N_DAM_MATCH(IB,4),K) - (EXFLUX-FXX-FYY) 469END IF 470END IF 471 #endif 472end do 473end do !!sigma loop 474 475 !--------------------------------------------------------------------------------- 476 ! Accumulate Fluxes at Boundary Nodes 477 !--------------------------------------------------------------------------------- 478 479 # if defined (MULTIPROCESSOR) 480if(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 !--------------------------------------------------------------------------------- 486do k=1,kbm1 487if(iobcn > 0) then 488do i=1,iobcn 489i1=i_obc_n(i) 490fflux_obc(i,k)=xflux_adv(i1,k) 491end do 492end if 493end do 494 495 !--------------------------------------------------------------------------------- 496 ! Calculate Vertical Advection Terms 497 !--------------------------------------------------------------------------------- 498 499do i=1,m 500 #if defined (WET_DRY) 501if(iswetn(i)*iswetnt(i) == 1) then 502 #endif 503 #if defined (THIN_DAM) 504if(IS_DAM(I)==1)then 505wvel(1:kb)=0.0_sp 506call calc_vflux(kbm1,f(i,1:kbm1),wvel(1:kb),vflux) 507else 508call calc_vflux(kbm1,f(i,1:kbm1),wts(i,1:kb),vflux) 509end if 510 #else 511call calc_vflux(kbm1,f(i,1:kbm1),wts(i,1:kb),vflux) 512 #endif 513 514do k=1,kbm1 515if(isonb(i) == 2) then 516xflux(i,k)= (vflux(k)-vflux(k+1))*art1(i)/dz(i,k) 517else 518xflux(i,k)=xflux(i,k)+ (vflux(k)-vflux(k+1))*art1(i)/dz(i,k) 519end if 520 #if defined (THIN_DAM) 521IF(IS_DAM(I)==1.AND.K<=KDAM(I))THEN 522tmpflx = (vflux(k)-vflux(k+1))*art1(i)/dz(i,k) 523IF(N_DAM_MATCH(I,1)==1)THEN 524XFLUX(N_DAM_MATCH(I,2),K) = XFLUX(N_DAM_MATCH(I,2),K)+tmpflx 525END IF 526IF(N_DAM_MATCH(I,1)==2)THEN 527XFLUX(N_DAM_MATCH(I,2),K) = XFLUX(N_DAM_MATCH(I,2),K)+tmpflx 528XFLUX(N_DAM_MATCH(I,3),K) = XFLUX(N_DAM_MATCH(I,3),K)+tmpflx 529END IF 530IF(N_DAM_MATCH(I,1)==3)THEN 531XFLUX(N_DAM_MATCH(I,2),K) = XFLUX(N_DAM_MATCH(I,2),K)+tmpflx 532XFLUX(N_DAM_MATCH(I,3),K) = XFLUX(N_DAM_MATCH(I,3),K)+tmpflx 533XFLUX(N_DAM_MATCH(I,4),K) = XFLUX(N_DAM_MATCH(I,4),K)+tmpflx 534END IF 535END IF 536 #endif 537end do 538 #if defined (WET_DRY) 539end if 540 #endif 541end do 542 543 !------------------------------------------------------- 544 !Point Source 545 !------------------------------------------------------- 546if(source)then!!user specified 547 548if(RIVER_TS_SETTING == 'calculated') then 549if(RIVER_INFLOW_LOCATION == 'node') then 550do j=1,numqbc 551jj=inodeq(j) 552fpoint=fdis(j) 553do k=1,kbm1 554xflux(jj,k)=xflux(jj,k) - qdis(j)*vqdist(j,k)*fpoint !/dz(jj,k) 555end do 556end do 557else if(RIVER_INFLOW_LOCATION == 'edge') then 558write(*,*)'scalar advection not setup for "edge" point source' 559stop 560end if 561end if 562 563else 564 565if(RIVER_TS_SETTING == 'calculated') then 566if(RIVER_INFLOW_LOCATION == 'node') then 567do j=1,numqbc 568jj=inodeq(j) 569do k=1,kbm1 570fpoint = f(jj,k) 571xflux(jj,k)=xflux(jj,k) - qdis(j)*vqdist(j,k)*fpoint !/dz(jj,k) 572end do 573end do 574else if(RIVER_INFLOW_LOCATION == 'edge') then 575write(*,*)'scalar advection not setup for "edge" point source' 576stop 577end if 578end if 579 580endif 581 !------------------------------------------------------------------------ 582 !Update Scalar Quantity 583 !------------------------------------------------------------------------ 584 585do i=1,m 586 #if defined (WET_DRY) 587if(iswetn(i)*iswetnt(i) == 1 )then 588 #endif 589do k=1,kbm1 590 #if !defined (THIN_DAM) 591fn(i,k)=(f(i,k)-xflux(i,k)/art1(i)*(deltat/dt(i)))*(dt(i)/dtfa(i)) 592 #else 593IF(IS_DAM(I)==1.AND.K<=KDAM(I))THEN 594fn(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)) 596ELSE 597fn(i,k)=(f(i,k)-xflux(i,k)/art1(i)*(deltat/dt(i)))*(dt(i)/dtfa(i)) 598END IF 599 #endif 600end do 601 #if defined (WET_DRY) 602else 603do k=1,kbm1 604fn(i,k)=f(i,k) 605end do 606end if 607 #endif 608end do 609 610return 611End Subroutine Adv_Scal 612 !==============================================================================| 613 614 !==============================================================================| 615 ! Vertical Diffusion of Scalar| 616 !==============================================================================| 617Subroutine Vdif_Scal(f,deltat) 618 619use mTridiagonal 620use all_vars 621 # if defined (THIN_DAM) 622use 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 627Implicit None 628Real(sp), intent(inout) :: f(0:mt,kb) 629Real(sp), intent(in) :: deltat 630!--local-------------------- 631integer:: i,k,ll 632real(sp) :: dsqrd,dfdz,visb 633real(sp) :: fsol(0:kb) 634 635 # if defined (THIN_DAM) 636real(sp) :: ftmp,stmp 637 # endif 638 639call init_tridiagonal(kb) 640 641Do i=1,m 642dsqrd = d(i)*d(i) 643 644!---------------------------------------------------------------- 645!Set up Diagonals of Matrix (lower=au,diag=bu,upper=cu) 646!---------------------------------------------------------------- 647 648 649!Surface 650au(1) = 0.0 651cu(1)=- deltat*(kh(i,2)+umol)/(dzz(i,1)*dz(i,1)*dsqrd) 652bu(1)=1.0 - cu(1) 653 654!Interior 655do k=2,kbm1-1 656au(k) =- deltat*(kh(i,k)+umol)/(dzz(i,k-1)*dz(i,k)*dsqrd) 657cu(k) =- deltat*(kh(i,k+1)+umol)/(dzz(i,k)*dz(i,k)*dsqrd) 658bu(k) = 1.0 - cu(k) - au(k) 659end do 660 661!Bottom 662au(kbm1) =- deltat*(kh(i,kbm1)+umol)/(dzz(i,kbm1-1)*dz(i,kbm1)*dsqrd) 663cu(kbm1) = 0.0 664bu(kbm1) = 1.0 - au(kbm1) 665 666!---------------------------------------------------------------- 667! Set up RHS forcing vector and boundary conditions 668!---------------------------------------------------------------- 669do k=1,kbm1 670du(k) = f(i,k) 671end do 672 673!Free Surface: No flux 674 675!Bottom: No flux 676 677 678!---------------------------------------------------------------- 679! Solve 680!---------------------------------------------------------------- 681 682call tridiagonal(kb,1,kbm1,fsol) 683 684!Transfer 685f(i,1:kbm1) = fsol(1:kbm1) 686 687End Do 688 689 #if defined (THIN_DAM) 690DO K=1,KBM1 691DO I=1,NODE_DAM1_N 692IF(K<=KDAM(I_NODE_DAM1_N(I,1)).AND.K<=KDAM(I_NODE_DAM1_N(I,2)) )THEN 693FTMP=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)) 695STMP=ART1(I_NODE_DAM1_N(I,1))+ART1(I_NODE_DAM1_N(I,2)) 696F(I_NODE_DAM1_N(I,1),K)=FTMP/STMP 697F(I_NODE_DAM1_N(I,2),K)=FTMP/STMP 698END IF 699END DO 700 701DO I=1,NODE_DAM2_N 702IF(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 704FTMP= 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)) 707STMP=ART1(I_NODE_DAM2_N(I,1))+ART1(I_NODE_DAM2_N(I,2)) & 708&+ART1(I_NODE_DAM2_N(I,3)) 709F(I_NODE_DAM2_N(I,1),K)=FTMP/STMP 710F(I_NODE_DAM2_N(I,2),K)=FTMP/STMP 711F(I_NODE_DAM2_N(I,3),K)=FTMP/STMP 712END IF 713END DO 714 715DO I=1,NODE_DAM3_N 716IF(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 718FTMP =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)) 722STMP =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)) 724F(I_NODE_DAM3_N(I,1),K)=FTMP/STMP 725F(I_NODE_DAM3_N(I,2),K)=FTMP/STMP 726F(I_NODE_DAM3_N(I,3),K)=FTMP/STMP 727F(I_NODE_DAM3_N(I,4),K)=FTMP/STMP 728END IF 729END DO 730END DO 731 #endif 732 733 734End Subroutine Vdif_Scal 735 736 737 !==============================================================================| 738 ! Set Point Source Conditions for Scalar Function| 739 !==============================================================================| 740 741Subroutine Bcond_Scal_PTsource(f,fn,fdis) 742 743 !------------------------------------------------------------------------------| 744use all_vars 745use bcs 746use mod_obcs 747implicit none 748real(sp), intent(in ), dimension(0:mt,kb):: f 749real(sp), intent(out), dimension(0:mt,kb):: fn 750real(sp), intent(in ), dimension(numqbc ):: fdis 751 !--local------------------------------------------- 752integer:: i,j,k,j1,j11,j22 753 !------------------------------------------------------------------------------| 754 755 756 !-------------------------------------------- 757 ! Set Source Terms 758 !-------------------------------------------- 759if(RIVER_TS_SETTING == 'specified') then 760if(numqbc > 0) then 761if(RIVER_INFLOW_LOCATION == 'node') then 762do i=1,numqbc 763j11=inodeq(i) 764do k=1,kbm1 765fn(j11,k)=fdis(i) 766end do 767end do 768else if(RIVER_INFLOW_LOCATION == 'edge') then 769do i=1,numqbc 770j11=n_icellq(i,1) 771j22=n_icellq(i,2) 772do k=1,kbm1 773fn(j11,k)=fdis(i) 774fn(j22,k)=fdis(i) 775end do 776end do 777end if 778end if 779end if 780 781return 782End Subroutine Bcond_Scal_PTSource 783 !==============================================================================| 784 !==============================================================================| 785 786 !==============================================================================| 787 !Set Boundary Conditions for Scalar Function on Open Boundary| 788 !==============================================================================| 789 790Subroutine Bcond_Scal_OBC(f,fn,fflux_obc,f_obc,deltat,alpha_nudge) 791 792 !------------------------------------------------------------------------------| 793use all_vars 794use bcs 795use mod_obcs 796implicit none 797real(sp), intent(in), dimension(0:mt,kb):: f 798real(sp), intent(inout), dimension(0:mt,kb):: fn 799real(sp), intent(in), dimension(iobcn+1,kbm1) :: fflux_obc 800real(sp), intent(in), dimension(iobcn) :: f_obc 801real(sp), intent(in):: deltat 802real(sp), intent(in):: alpha_nudge 803 !--local------------------------------------------- 804real(sp) :: f2d,f2d_next,f2d_obc,xflux2d,tmp 805integer:: i,j,k,j1,j11,j22 806 !------------------------------------------------------------------------------| 807 808 !-------------------------------------------- 809 ! Set Scalar Value on Open Boundary 810 !-------------------------------------------- 811if(iobcn > 0) then 812do i=1,iobcn 813j=i_obc_n(i) 814j1=next_obc(i) 815f2d=0.0_sp 816f2d_next=0.0_sp 817xflux2d=0.0_sp 818do k=1,kbm1 819f2d=f2d+f(j,k)*dz(j,k) 820f2d_next=f2d_next+fn(j1,k)*dz(j1,k) 821xflux2d=xflux2d+fflux_obc(i,k)*dz(j,k) 822end do 823 824if(uard_obcn(i) > 0.0_sp) then 825tmp=xflux2d+f2d*uard_obcn(i) 826f2d_obc=(f2d*dt(j)-tmp*deltat/art1(j))/d(j) 827do k=1,kbm1 828fn(j,k)=fn(j1,k) !f2d_obc+(fn(j1,k)-f2d_next) 829end do 830else 831do k=1,kbm1 832fn(j,k) = f(j,k)-alpha_nudge*(f(j,k)-f_obc(i)) 833end do 834end if 835end do 836endif 837 838return 839End Subroutine Bcond_Scal_OBC 840 !==============================================================================| 841 !==============================================================================| 842 843Subroutine fct_sed(f,fn) 844!==============================================================================| 845USE ALL_VARS 846USE MOD_UTILS 847USE BCS 848USE MOD_OBCS 849IMPLICIT NONE 850real(sp), intent(inout), dimension(0:mt,kb):: fn 851real(sp), intent(in), dimension(0:mt,kb):: f 852REAL(SP):: SMAX,SMIN 853INTEGER :: I,J,K,K1 854!==============================================================================| 855IF(DBG_SET(DBG_SBR)) WRITE(IPT,*)"Start: fct_sed" 856 857nodes: DO I=1,M 858 859! SKIP OPEN BOUNDARY NODES 860IF(IOBCN > 0)THEN 861DO J=1,IOBCN 862IF(I == I_OBC_N(J)) CYCLE nodes 863END DO 864END IF 865 866! SKIP RIVER INFLOW POINTS 867IF(NUMQBC > 0)THEN 868DO J=1,NUMQBC 869IF(RIVER_INFLOW_LOCATION == 'node')THEN 870IF(I == INODEQ(J)) CYCLE nodes 871END IF 872IF(RIVER_INFLOW_LOCATION == 'edge')THEN 873IF(I == N_ICELLQ(J,1) .OR. I == N_ICELLQ(J,2)) CYCLE nodes 874END IF 875END DO 876END IF 877 878! SKIP GROUND WATER INFLOW POINTS 879IF(BFWDIS(I) .GT. 0.0_SP .and. GROUNDWATER_SALT_ON) CYCLE nodes 880 881K1 = 1 882IF(PRECIPITATION_ON) K1 = 2 883 !DO K=1,KBM1 884DO K=K1,KBM1 885SMAX = MAXVAL(f(NBSN(I,1:NTSN(I)),K)) 886SMIN = MINVAL(f(NBSN(I,1:NTSN(I)),K)) 887 888IF(K == 1)THEN 889SMAX = MAX(SMAX,(f(I,K)*DZ(I,K+1)+f(I,K+1)*DZ(I,K))/& 890(DZ(I,K)+DZ(I,K+1))) 891SMIN = MIN(SMIN,(f(I,K)*DZ(I,K+1)+f(I,K+1)*DZ(I,K))/& 892(DZ(I,K)+DZ(I,K+1))) 893ELSE IF(K == KBM1)THEN 894SMAX = MAX(SMAX,(f(I,K)*DZ(I,K-1)+f(I,K-1)*DZ(I,K))/& 895(DZ(I,K)+DZ(I,K-1))) 896SMIN = MIN(SMIN,(f(I,K)*DZ(I,K-1)+f(I,K-1)*DZ(I,K))/& 897(DZ(I,K)+DZ(I,K-1))) 898ELSE 899SMAX = 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))) 903SMIN = 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))) 907END IF 908 909IF(SMIN-fn(I,K) > 0.0_SP)fn(I,K) = SMIN 910IF(fn(I,K)-SMAX > 0.0_SP)fn(I,K) = SMAX 911 912END DO 913END DO nodes 914 915WHERE(fn < 0.0_SP)fn=0.0_SP 916 917IF(DBG_SET(DBG_SBR)) WRITE(IPT,*)"End: fct_sed" 918End 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 !========================================================================== 930Subroutine Calc_VFlux(n,c,w,flux) 931use mod_prec 932implicit none 933integer , intent(in ) :: n 934real(sp), intent(in ) ::c(n) 935real(sp), intent(in ) ::w(n+1) 936real(sp), intent(out) ::flux(n+1) 937real(sp) :: conv(n+1),diss(n+1) 938real(sp) :: cin(-1:n+2) 939real(sp) :: dis4 940integer:: i 941 942!transfer to working array 943cin(1:n) = c(1:n) 944 945!surface bcs (no flux) 946cin(0)=-cin(1) 947cin(-1) =-cin(2) 948 949!bottom bcs (no flux) 950cin(n+1) = -cin(n) 951cin(n+2) = -cin(n-1) 952 953!flux computation 954do i=1,n+1 955dis4= .5*abs(w(i)) 956conv(i) = w(i)*(cin(i)+cin(i-1))/2. 957diss(i) = dis4*(cin(i)-cin(i-1)-lim(cin(i+1)-cin(i),cin(i-1)-cin(i-2))) 958flux(i) = conv(i)+diss(i) 959end do 960 961End Subroutine Calc_VFlux 962 963 !========================================================================== 964 ! Calculate LED Limiter L(u,v) 965 !========================================================================== 966Function Lim(a,b) 967use mod_prec 968real(sp) lim,a,b 969real(sp) q,R 970real(sp) eps 971eps = epsilon(eps) 972 973q = 0. !1st order 974q = 1. !minmod 975q = 2. !van leer 976 977R = abs((a-b)/(abs(a)+abs(b)+eps) )**q 978lim = .5*(1-R)*(a+b) 979 980End Function Lim 981 982 983 End Module Scalar

View Code
为了更快速的了解计算过程，自己设置了一个只有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))*viscoffxx=-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+fyyxflux(ia,k)=xflux(ia,k)+exflux xflux(ib,k)=xflux(ib,k)-exfluxfn(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的二阶精度计算，只需按照泰勒公式计算即可