Actual source code: ex5.c
1: static char help[] = "Nonlinear, time-dependent. Developed from radiative_surface_balance.c - work with petsc-dev \n";
2: /*
3: Contributed by Steve Froehlich, Illinois Institute of Technology
5: Usage:
6: mpiexec -n <np> ./ex5 [options]
7: ./ex5 -help [view petsc options]
8: ./ex5 -ts_type sundials -ts_sundials_monitor_steps -pc_type lu -ts_view
9: ./ex5 -da_grid_x 20 -da_grid_y 20 -log_summary
10: ./ex5 -drawcontours -draw_pause 0.1 -draw_fields 0,1,2,3,4
11: */
13: /*
14: -----------------------------------------------------------------------
16: Governing equations:
18: R = s*(Ea*Ta^4 - Es*Ts^4)
19: SH = p*Cp*Ch*wind*(Ta - Ts)
20: LH = p*L*Ch*wind*B(q(Ta) - q(Ts))
21: G = k*(Tgnd - Ts)/dz
23: Fnet = R + SH + LH + G
25: du/dt = -u*(du/dx) - v*(du/dy) - 2*omeg*sin(lat)*v - (1/p)*(dP/dx)
26: dv/dt = -u*(dv/dx) - v*(dv/dy) + 2*omeg*sin(lat)*u - (1/p)*(dP/dy)
27: dTs/dt = Fnet/(Cp*dz) - Div([u*Ts, v*Ts]) + D*Lap(Ts)
28: = Fnet/(Cs*dz) - u*(dTs/dx) - v*(dTs/dy) + D*(Ts_xx + Ts_yy)
29: dp/dt = -Div([u*p,v*p])
30: = - u*dp/dx - v*dp/dy
31: dTa/dt = Fnet/Cp
33: Equation of State:
35: P = p*R*Ts
37: -----------------------------------------------------------------------
39: Program considers the evolution of a two dimensional atmosphere from
40: sunset to sunrise. There are two components:
41: 1. Surface energy balance model to compute diabatic dT (Fnet)
42: 2. Dynamical model using simplified primitive equations
44: Program is to be initiated at sunset and run to sunrise.
46: Inputs are:
47: Surface temperature
48: Dew point temperature
49: Air temperature
50: Temperature at cloud base (if clouds are present)
51: Fraction of sky covered by clouds
52: Wind speed
53: Precipitable water in centimeters
54: Wind direction
56: Inputs are are read in from the text file ex5_control.txt. To change an
57: input value use ex5_control.txt.
59: Solvers:
60: Backward Euler = default solver
61: Sundials = fastest and most accurate, requires Sundials libraries
63: This model is under development and should be used only as an example
64: and not as a predictive weather model.
65: */
67: #include <petscts.h>
68: #include <petscdmda.h>
70: #define SIG 0.000000056703 //stefan-boltzmann constant
71: #define EMMSFC 1 //absorption-emission constant for surface
72: #define TIMESTEP 1 //amount of time(seconds) that passes before new flux is calculated
74: /* variables of interest to be solved at each grid point */
75: typedef struct {
76: PetscScalar Ts,Ta; // surface and air temperature
77: PetscScalar u,v; // wind speed
78: PetscScalar p; // density
79: } Field;
81: /* User defined variables. Used in solving for variables of interest */
82: typedef struct {
83: DM da; //grid
84: PetscScalar csoil; //heat constant for layer
85: PetscScalar dzlay; //thickness of top soil layer
86: PetscScalar emma; //emission parameter
87: PetscScalar wind; //wind speed
88: PetscScalar dewtemp; //dew point temperature (moisture in air)
89: PetscScalar pressure1; //sea level pressure
90: PetscScalar airtemp; //temperature of air near boundary layer inversion
91: PetscScalar Ts; //temperature at the surface
92: PetscScalar fract; //fraction of sky covered by clouds
93: PetscScalar Tc; //temperature at base of lowest cloud layer
94: PetscScalar lat; //Latitude in degrees
95: PetscScalar init; //initialization scenario
96: PetscScalar deep_grnd_temp;//temperature of ground under top soil surface layer
97: } AppCtx;
99: /* Struct for visualization */
100: typedef struct {
101: PetscBool drawcontours; /* flag - 1 indicates drawing contours */
102: PetscViewer drawviewer;
103: } MonitorCtx;
106: /* Inputs read in from text file */
107: struct in {
108: PetscScalar Ts; //surface temperature
109: PetscScalar Td; //dewpoint temperature
110: PetscScalar Tc; //temperature of cloud base
111: PetscScalar fr; //fraction of sky covered by clouds
112: PetscScalar wnd; //wind speed
113: PetscScalar Ta; //air temperature
114: PetscScalar pwt; //precipitable water
115: PetscScalar wndDir; //wind direction
116: PetscScalar lat; //latitude
117: PetscReal time; //time in hours
118: PetscScalar init;
119: };
121: //functions
143: int main(int argc,char **argv)
144: {
145: PetscErrorCode ierr;
146: int time; //amount of loops
147: struct in put;
148: PetscScalar rh; //relative humidity
149: PetscScalar x; //memory varialbe for relative humidity calculation
150: PetscScalar deep_grnd_temp; //temperature of ground under top soil surface layer
152: PetscScalar emma; //absorption-emission constant for air
153: PetscScalar pressure1 = 101300; //surface pressure
154: PetscScalar mixratio; //mixing ratio
155: PetscScalar airtemp; //temperature of air near boundary layer inversion
156: PetscScalar dewtemp; //dew point temperature
157: PetscScalar sfctemp; //temperature at surface
158: PetscScalar pwat; //total column precipitable water
159: PetscScalar cloudTemp; //temperature at base of cloud
160: AppCtx user; /* user-defined work context */
161: MonitorCtx usermonitor; /* user-defined monitor context */
162: PetscMPIInt rank,size;
164: PetscInitialize(&argc,&argv,(char *)0,help);
165: MPI_Comm_size(PETSC_COMM_WORLD,&size);
166: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
168: /* Inputs */
169: readinput(&put);
171: sfctemp = put.Ts;
172: dewtemp = put.Td;
173: cloudTemp = put.Tc;
174: airtemp = put.Ta;
175: pwat = put.pwt;
177: if (!rank) PetscPrintf(PETSC_COMM_SELF,"Initial Temperature = %g\n",sfctemp); //input surface temperature
179: deep_grnd_temp = sfctemp - 10; //set underlying ground layer temperature
180: emma = emission(pwat); //accounts for radiative effects of water vapor
182: /* Converts from Fahrenheit to Celsuis */
183: sfctemp = fahr_to_cel(sfctemp);
184: airtemp = fahr_to_cel(airtemp);
185: dewtemp = fahr_to_cel(dewtemp);
186: cloudTemp = fahr_to_cel(cloudTemp);
187: deep_grnd_temp = fahr_to_cel(deep_grnd_temp);
189: /* Converts from Celsius to Kelvin */
190: sfctemp +=273;
191: airtemp +=273;
192: dewtemp +=273;
193: cloudTemp +=273;
194: deep_grnd_temp +=273;
196: /* Calculates initial relative humidity */
197: x = calcmixingr(dewtemp,pressure1);
198: mixratio = calcmixingr(sfctemp,pressure1);
199: rh = (x/mixratio)*100;
201: if (!rank){printf("Initial RH = %.1f percent\n\n",rh);} //prints initial relative humidity
203: time = 3600*put.time; //sets amount of timesteps to run model
205: /* Use PETSc TS solver */
206: /*------------------------------------------*/
207: TS ts;
208: DM da;
209: Vec T,rhs; /* solution vector */
210: Mat J; /* Jacobian matrix */
211: PetscReal ftime,dt;
212: PetscInt steps,dof=5;
214: /*Create grid*/
215: DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_PERIODIC,DMDA_BOUNDARY_PERIODIC,DMDA_STENCIL_STAR,-20,-20,
216: PETSC_DECIDE,PETSC_DECIDE,dof,1,PETSC_NULL,PETSC_NULL,&da);
218: /*Define output window for each variable of interest*/
219: DMDASetFieldName(da,0,"Ts");
220: DMDASetFieldName(da,1,"Ta");
221: DMDASetFieldName(da,2,"u");
222: DMDASetFieldName(da,3,"v");
223: DMDASetFieldName(da,4,"p");
225: /*set values for appctx*/
226: user.da = da;
227: user.Ts = sfctemp;
228: user.fract = put.fr; //fraction of sky covered by clouds
229: user.dewtemp = dewtemp; //dew point temperature (mositure in air)
230: user.csoil = 2000000; //heat constant for layer
231: user.dzlay = 0.08; //thickness of top soil layer
232: user.emma = emma; //emission parameter
233: user.wind = put.wnd; //wind spped
234: user.pressure1 = pressure1; //sea level pressure
235: user.airtemp = airtemp; //temperature of air near boundar layer inversion
236: user.Tc = cloudTemp; //temperature at base of lowest cloud layer
237: user.init = put.init; //user chosen initiation scenario
238: user.lat = 70*0.0174532; //converts latitude degrees to latitude in radians
239: user.deep_grnd_temp = deep_grnd_temp; //temp in lowest ground layer
241: /*set values for MonitorCtx*/
242: usermonitor.drawcontours = PETSC_FALSE;
243: PetscOptionsHasName(PETSC_NULL,"-drawcontours",&usermonitor.drawcontours);
244: if (usermonitor.drawcontours){
245: PetscReal bounds[] = {1000.0,-1000., -1000.,-1000., 1000.,-1000., 1000.,-1000., 1000,-1000, 100700,100800};
246: PetscViewerDrawOpen(PETSC_COMM_WORLD,0,0,0,0,300,300,&usermonitor.drawviewer);
247: PetscViewerDrawSetBounds(usermonitor.drawviewer,dof,bounds);
248: }
250: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
251: Extract global vectors from DA;
252: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
253: DMCreateGlobalVector(da,&T);
254: VecDuplicate(T,&rhs); //r: vector to put the computed right hand side
256: TSCreate(PETSC_COMM_WORLD,&ts);
257: TSSetProblemType(ts,TS_NONLINEAR);
258: TSSetType(ts,TSBEULER);
259: TSSetRHSFunction(ts,rhs,RhsFunc,&user);
261: /* Set Jacobian evaluation routine - use coloring to compute finite difference Jacobian efficiently */
262: PetscBool use_coloring=PETSC_TRUE;
263: MatFDColoring matfdcoloring=0;
264: DMGetMatrix(da,MATAIJ,&J);
265: if (use_coloring){
266: ISColoring iscoloring;
267: DMGetColoring(da,IS_COLORING_GLOBAL,MATAIJ,&iscoloring);
268: MatFDColoringCreate(J,iscoloring,&matfdcoloring);
269: MatFDColoringSetFromOptions(matfdcoloring);
270: ISColoringDestroy(&iscoloring);
271: MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))RhsFunc,&user);
272: TSSetRHSJacobian(ts,J,J,TSDefaultComputeJacobianColor,matfdcoloring);
273: } else {
274: TSSetRHSJacobian(ts,J,J,TSDefaultComputeJacobian,PETSC_NULL);
275: }
277: /*Define what to print for ts_monitor option*/
278: TSMonitorSet(ts,Monitor,&usermonitor,PETSC_NULL);
279: FormInitialSolution(da,T,&user);
280: dt = TIMESTEP; /* initial time step */
281: ftime = TIMESTEP*time;
282: if (!rank){printf("time %d, ftime %g hour, TIMESTEP %g\n",time,ftime/3600,dt);}
283: TSSetInitialTimeStep(ts,0.0,dt);
284: TSSetDuration(ts,time,ftime);
285: TSSetSolution(ts,T);
287: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
288: Set runtime options
289: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
290: TSSetFromOptions(ts);
292: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
293: Solve nonlinear system
294: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
295: TSSolve(ts,T,&ftime);
296: TSGetTimeStepNumber(ts,&steps);
297: if (!rank){PetscPrintf(PETSC_COMM_WORLD,"Solution T after %g hours %d steps\n",ftime/3600,steps);}
300: if (matfdcoloring){MatFDColoringDestroy(&matfdcoloring);}
301: if (usermonitor.drawcontours){
302: PetscViewerDestroy(&usermonitor.drawviewer);
303: }
304: MatDestroy(&J);
305: VecDestroy(&T);
306: VecDestroy(&rhs);
307: TSDestroy(&ts);
308: DMDestroy(&da);
310: PetscFinalize();
311: return 0;
312: }
313: /*****************************end main program********************************/
314: /*****************************************************************************/
315: /*****************************************************************************/
316: /*****************************************************************************/
319: PetscErrorCode calcfluxs(PetscScalar sfctemp, PetscScalar airtemp, PetscScalar emma, PetscScalar fract, PetscScalar cloudTemp, PetscScalar* flux)
320: {
322: *flux = SIG*((EMMSFC*emma*pow(airtemp,4)) + (EMMSFC*fract*(1 - emma)*pow(cloudTemp,4)) - (EMMSFC*pow(sfctemp,4))); //calculates flux using Stefan-Boltzmann relation
323: return(0);
324: }
328: PetscErrorCode calcfluxa(PetscScalar sfctemp, PetscScalar airtemp, PetscScalar emma, PetscScalar* flux) //this function is not currently called upon
329: {
330: PetscScalar emm = 0.001;
332: *flux = SIG*(- emm*(pow(airtemp,4))); //calculates flux usinge Stefan-Boltzmann relation
333: return(0);
334: }
337: PetscErrorCode sensibleflux(PetscScalar sfctemp, PetscScalar airtemp, PetscScalar wind, PetscScalar* sheat)
338: {
339: PetscScalar density = 1; //air density
340: PetscScalar Cp = 1005; //heat capicity for dry air
341: PetscScalar wndmix; //temperature change from wind mixing: wind*Ch
345: wndmix = 0.0025 + 0.0042*wind; //regression equation valid for neutral and stable BL
346: *sheat = density*Cp*wndmix*(airtemp - sfctemp); //calculates sensible heat flux
348: return(0);
349: }
353: PetscErrorCode latentflux(PetscScalar sfctemp, PetscScalar dewtemp, PetscScalar wind, PetscScalar pressure1, PetscScalar* latentheat)
354: {
355: PetscScalar density = 1; //density of dry air
356: PetscScalar q; //actual specific humitity
357: PetscScalar qs; //saturation specific humidity
358: PetscScalar wndmix; //temperature change from wind mixing: wind*Ch
359: PetscScalar beta = .4; //moisture availability
360: PetscScalar mr ; //mixing ratio
361: PetscScalar lhcnst; //latent heat of vaporization constant = 2501000 J/kg at 0c
362: //latent heat of saturation const = 2834000 J/kg
363: //latent heat of fusion const = 333700 J/kg
366: wind = mph2mpers(wind); //converts wind from mph to meters per second
367: wndmix = 0.0025 + 0.0042*wind; //regression equation valid for neutral BL
368: lhcnst = Lconst(sfctemp); //calculates latent heat of evaporation
369: mr = calcmixingr(sfctemp,pressure1);//calculates saturation mixing ratio
370: qs = calc_q(mr); //calculates saturation specific humidty
371: mr = calcmixingr(dewtemp,pressure1); //calculates mixing ratio
372: q = calc_q(mr); //calculates specific humidty
374: *latentheat = density*wndmix*beta*lhcnst*(q - qs); //calculates latent heat flux
375: return(0);
376: }
380: PetscErrorCode potential_temperature(PetscScalar temp, PetscScalar pressure1, PetscScalar pressure2, PetscScalar sfctemp, PetscScalar* pottemp)
381: {
382: PetscScalar kdry; //poisson constant for dry atmosphere
383: PetscScalar kmoist; //poisson constant for moist atmosphere
384: PetscScalar pavg; //average atmospheric pressure
385: PetscScalar mixratio;//mixing ratio
388: mixratio = calcmixingr(sfctemp,pressure1);
390: /*initialize poisson constant */
391: kdry = 0.2854;
392: kmoist = 0.2854*(1 - 0.24*mixratio);
394: pavg = ((0.7*pressure1)+pressure2)/2; //calculates simple average press
395: *pottemp = temp*(pow((pressure1/pavg),kdry)); //calculates potential temperature
397: return(0);
398: }
400: {
401: PetscScalar e; //vapor pressure
402: PetscScalar mixratio; //mixing ratio
404: dtemp = dtemp - 273; //converts from Kelvin to Celsuis
405: e = 6.11*(pow(10,((7.5*dtemp)/(237.7+dtemp)))); //converts from dew point temp to vapor pressure
406: e = e*100; //converts from hPa to Pa
407: mixratio = (0.622*e)/(pressure1 - e); //computes mixing ratio
408: mixratio = mixratio*1; //convert to g/Kg
410: return mixratio;
411: }
413: {
414: PetscScalar specific_humidity; //define specific humidity variable
415: specific_humidity = rv/(1 + rv); //calculates specific humidity
416: return specific_humidity;
417: }
421: PetscErrorCode calc_gflux(PetscScalar sfctemp, PetscScalar deep_grnd_temp, PetscScalar* Gflux)
422: {
423: PetscScalar k; //thermal conductivity parameter
424: PetscScalar n = 0.38; //value of soil porosity
425: PetscScalar dz = 1; //depth of layer between soil surface and deep soil layer
426: PetscScalar unit_soil_weight = 2700; //unit soil weight in kg/m^3
430: k = ((0.135*(1-n)*unit_soil_weight) + 64.7)/(unit_soil_weight - (0.947*(1-n)*unit_soil_weight)); //dry soil conductivity
431: *Gflux = (k*(deep_grnd_temp - sfctemp)/dz); //calculates flux from deep ground layer
432: return(0);
433: }
437: {
438: PetscScalar emma;
440: emma = 0.725 + 0.17*log10(pwat);
442: return emma;
443: }
445: {
446: PetscScalar emma = 0;
448: /*modifies radiative balance depending on cloud cover */
449: if (fract >= 0.9)
450: emma = 1;
451: else if (0.9 > fract && fract >= 0.8)
452: emma = 0.9;
453: else if (0.8 > fract && fract >= 0.7)
454: emma = 0.85;
455: else if (0.7 > fract && fract >= 0.6)
456: emma = 0.75;
457: else if (0.6 > fract && fract >= 0.5)
458: emma = 0.65;
459: else if (0.4 > fract && fract >= 0.3)
460: emma = emma*1.086956;
461: return emma;
462: }
464: {
465: PetscScalar Lheat;
466: sfctemp -=273; //converts from kelvin to celsius
467: Lheat = 4186.8*(597.31 - 0.5625*sfctemp); //calculates latent heat constant
468: return Lheat;
469: }
471: {
472: wind = ((wind*1.6*1000)/3600); //converts wind from mph to meters per second
473: return wind;
474: }
476: {
477: temp = (5*(temp-32))/9; //converts from farhrenheit to celsuis
478: return temp;
479: }
481: {
482: temp = ((temp*9)/5) + 32; //converts from celsuis to farhrenheit
483: return temp;
484: }
485: void readinput(struct in *put)
486: {
487: int i;
488: char x;
489: FILE *ifp;
491: ifp = fopen("ex5_control.txt", "r");
493: for (i=0;i<110;i++)
494: fscanf(ifp, "%c", &x);
495: fscanf(ifp, "%lf", &put->Ts);
497: for (i=0;i<43;i++)
498: fscanf(ifp, "%c", &x);
499: fscanf(ifp, "%lf", &put->Td);
501: for (i=0;i<43;i++)
502: fscanf(ifp, "%c", &x);
503: fscanf(ifp, "%lf", &put->Ta);
505: for (i=0;i<43;i++)
506: fscanf(ifp, "%c", &x);
507: fscanf(ifp, "%lf", &put->Tc);
509: for (i=0;i<43;i++)
510: fscanf(ifp, "%c", &x);
511: fscanf(ifp, "%lf", &put->fr);
513: for (i=0;i<43;i++)
514: fscanf(ifp, "%c", &x);
515: fscanf(ifp, "%lf", &put->wnd);
517: for (i=0;i<43;i++)
518: fscanf(ifp, "%c", &x);
519: fscanf(ifp, "%lf", &put->pwt);
521: for (i=0;i<43;i++)
522: fscanf(ifp, "%c", &x);
523: fscanf(ifp, "%lf", &put->wndDir);
525: for (i=0;i<43;i++)
526: fscanf(ifp, "%c", &x);
527: fscanf(ifp, "%lf", &put->time);
529: for (i=0;i<63;i++)
530: fscanf(ifp, "%c", &x);
531: fscanf(ifp, "%lf", &put->init);
532: }
534: /* ------------------------------------------------------------------- */
537: PetscErrorCode FormInitialSolution(DM da,Vec Xglobal,void *ctx)
538: {
540: AppCtx *user = (AppCtx*)ctx; /* user-defined application context */
541: PetscInt i,j,xs,ys,xm,ym,Mx,My;
542: Field **X;
543: PetscScalar deltT;
544: PetscReal hx,hy;
545: FILE *ifp;
546: FILE *ofp;
549: ofp = fopen("swing", "w");
550: ifp = fopen("grid.in", "r");
551: deltT = 0.8;
553: DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
554: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
556: hx = 1/(PetscReal)(Mx-1);
557: hy = 1/(PetscReal)(My-1);
559: /* Get pointers to vector data */
560: DMDAVecGetArray(da,Xglobal,&X);
562: /* Get local grid boundaries */
563: DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
565: /* Compute function over the locally owned part of the grid */
567: if (user->init == 1) {
568: for (j=ys; j<ys+ym; j++) {
569: for (i=xs; i<xs+xm; i++) {
570: X[j][i].Ts = user->Ts - i*0.0001;
571: X[j][i].Ta = X[j][i].Ts - 5;
572: X[j][i].u = 0;
573: X[j][i].v = 0;
574: X[j][i].p = 1.25;
575: if ((j == 5 || j == 6) && (i == 4 || i == 5)) X[j][i].p += 0.00001;
576: if ((j == 5 || j == 6) && (i == 12 || i == 13)) X[j][i].p += 0.00001;
577: }
578: }
579: }
581: else {
582: for (j=ys; j<ys+ym; j++) {
583: for (i=xs; i<xs+xm; i++) {
584: X[j][i].Ts = user->Ts;
585: X[j][i].Ta = X[j][i].Ts - 5;
586: X[j][i].u = 0;
587: X[j][i].v = 0;
588: X[j][i].p = 1.25;
589: }
590: }
591: }
593: /* Restore vectors */
594: DMDAVecRestoreArray(da,Xglobal,&X);
595: return(0);
596: }
600: /*
601: RhsFunc - Evaluates nonlinear function F(u).
603: Input Parameters:
604: . ts - the TS context
605: . t - current time
606: . Xglobal - input vector
607: . F - output vector
608: . ptr - optional user-defined context, as set by SNESSetFunction()
610: Output Parameter:
611: . F - rhs function vector
612: */
613: PetscErrorCode RhsFunc(TS ts,PetscReal t,Vec Xglobal,Vec F,void *ctx)
614: {
615: AppCtx *user = (AppCtx*)ctx; /* user-defined application context */
616: DM da = user->da;
618: PetscInt i,j,Mx,My,xs,ys,xm,ym;
619: PetscReal dhx,dhy;
620: Vec localT;
621: Field **X,**Frhs; //structures that contain variables of interest and left hand side of governing equations respectively
622: PetscScalar csoil = user->csoil; //heat constant for layer
623: PetscScalar dzlay = user->dzlay; //thickness of top soil layer
624: PetscScalar emma = user->emma; //emission parameter
625: PetscScalar wind = user->wind; //wind speed
626: PetscScalar dewtemp = user->dewtemp; //dew point temperature (moisture in air)
627: PetscScalar pressure1 = user->pressure1; //sea level pressure
628: PetscScalar airtemp = user->airtemp; //temperature of air near boundary layer inversion
629: PetscScalar fract = user->fract; //fraction of the sky covered by clouds
630: PetscScalar Tc = user->Tc; //temperature at base of lowest cloud layer
631: PetscScalar lat = user->lat; //latitude
632: PetscScalar Cp = 1005.7; //specific heat of air at constant pressure
633: PetscScalar Rd = 287.058; //gas constant for dry air
634: PetscScalar diffconst = 1000; //diffusion coefficient
635: PetscScalar f = 2*0.0000727*sin(lat); //coriolis force
636: PetscScalar deep_grnd_temp = user->deep_grnd_temp; //temp in lowest ground layer
637: PetscScalar Ts,u,v,p,P;
638: PetscScalar u_abs,u_plus,u_minus,v_abs,v_plus,v_minus;
640: PetscScalar sfctemp1,fsfc1,Ra;
641: PetscScalar sheat; //sensible heat flux
642: PetscScalar latentheat; //latent heat flux
643: PetscScalar groundflux; //flux from conduction of deep ground layer in contact with top soil
644: PetscInt xend,yend;
647: DMGetLocalVector(da,&localT);
648: DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
649: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
651: dhx = (PetscReal)(Mx-1)/(5000*(Mx-1)); // dhx = 1/dx; assume 2D space domain: [0.0, 1.e5] x [0.0, 1.e5]
652: dhy = (PetscReal)(My-1)/(5000*(Mx-1)); // dhy = 1/dy;
655: /*
656: Scatter ghost points to local vector,using the 2-step process
657: DAGlobalToLocalBegin(),DAGlobalToLocalEnd().
658: By placing code between these two statements, computations can be
659: done while messages are in transition.
660: */
661: DMGlobalToLocalBegin(da,Xglobal,INSERT_VALUES,localT);
662: DMGlobalToLocalEnd(da,Xglobal,INSERT_VALUES,localT);
664: /* Get pointers to vector data */
665: DMDAVecGetArray(da,localT,&X);
666: DMDAVecGetArray(da,F,&Frhs);
668: /* Get local grid boundaries */
669: DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
671: /* Compute function over the locally owned part of the grid */
672: /* the interior points */
673: xend=xs+xm; yend=ys+ym;
674: for (j=ys; j<yend; j++) {
675: for (i=xs; i<xend; i++) {
676: Ts = X[j][i].Ts; u = X[j][i].u; v = X[j][i].v; p = X[j][i].p; //P = X[j][i].P;
678: sfctemp1 = (double)Ts;
679: sfctemp1 = (double)X[j][i].Ts;
680: calcfluxs(sfctemp1,airtemp,emma,fract,Tc,&fsfc1); //calculates surface net radiative flux
681: sensibleflux(sfctemp1,airtemp,wind,&sheat); //calculate sensible heat flux
682: latentflux(sfctemp1,dewtemp,wind,pressure1,&latentheat);//calculates latent heat flux
683: calc_gflux(sfctemp1,deep_grnd_temp,&groundflux); //calculates flux from earth below surface soil layer by conduction
684: calcfluxa(sfctemp1,airtemp,emma,&Ra); //Calculates the change in downward radiative flux
685: fsfc1 = fsfc1 + latentheat + sheat + groundflux; //adds radiative, sensible heat, latent heat, and ground heat flux yielding net flux
687: /* convective coefficients for upwinding */
688: u_abs = PetscAbsScalar(u);
689: u_plus = .5*(u + u_abs); // u if u>0; 0 if u<0
690: u_minus = .5*(u - u_abs); // u if u <0; 0 if u>0
692: v_abs = PetscAbsScalar(v);
693: v_plus = .5*(v + v_abs); // v if v>0; 0 if v<0
694: v_minus = .5*(v - v_abs); // v if v <0; 0 if v>0
696: /* Solve governing equations */
697: P = p*Rd*Ts;
699: /* du/dt -> time change of east-west component of the wind */
700: Frhs[j][i].u = - u_plus*(u - X[j][i-1].u)*dhx - u_minus*(X[j][i+1].u - u)*dhx // - u(du/dx)
701: - v_plus*(u - X[j-1][i].u)*dhy - v_minus*(X[j+1][i].u - u)*dhy // - v(du/dy)
702: -(Rd/p)*(Ts*(X[j][i+1].p - X[j][i-1].p)*0.5*dhx + p*0*(X[j][i+1].Ts - X[j][i-1].Ts)*0.5*dhx) // -(R/p)[Ts(dp/dx)+ p(dTs/dx)]
703: // -(1/p)*(X[j][i+1].P - X[j][i-1].P)*dhx
704: + f*v;
706: /* dv/dt -> time change of north-south component of the wind */
707: Frhs[j][i].v = - u_plus*(v - X[j][i-1].v)*dhx - u_minus*(X[j][i+1].v - v)*dhx // - u(dv/dx)
708: - v_plus*(v - X[j-1][i].v)*dhy - v_minus*(X[j+1][i].v - v)*dhy // - v(dv/dy)
709: -(Rd/p)*(Ts*(X[j+1][i].p - X[j-1][i].p)*0.5*dhy + p*0*(X[j+1][i].Ts - X[j-1][i].Ts)*0.5*dhy) // -(R/p)[Ts(dp/dy)+ p(dTs/dy)]
710: // -(1/p)*(X[j+1][i].P - X[j-1][i].P)*dhy
711: -f*u;
713: /* dT/dt -> time change of temperature */
714: Frhs[j][i].Ts = (fsfc1/(csoil*dzlay)) // Fnet/(Cp*dz) diabatic change in T
715: -u_plus*(Ts - X[j][i-1].Ts)*dhx - u_minus*(X[j][i+1].Ts - Ts)*dhx // - u*(dTs/dx) advection x
716: -v_plus*(Ts - X[j-1][i].Ts)*dhy - v_minus*(X[j+1][i].Ts - Ts)*dhy // - v*(dTs/dy) advection y
717: + diffconst*((X[j][i+1].Ts - 2*Ts + X[j][i-1].Ts)*dhx*dhx // + D(Ts_xx + Ts_yy) diffusion
718: + (X[j+1][i].Ts - 2*Ts + X[j-1][i].Ts)*dhy*dhy);
720: /* dp/dt -> time change of */
721: Frhs[j][i].p = -u_plus*(p - X[j][i-1].p)*dhx - u_minus*(X[j][i+1].p - p)*dhx // - u*(dp/dx)
722: -v_plus*(p - X[j-1][i].p)*dhy - v_minus*(X[j+1][i].p - p)*dhy; // - v*(dp/dy)
724: Frhs[j][i].Ta = Ra/Cp; //dTa/dt time change of air temperature
725: }
726: }
728: /* Restore vectors */
729: DMDAVecRestoreArray(da,localT,&X);
730: DMDAVecRestoreArray(da,F,&Frhs);
731: DMRestoreLocalVector(da,&localT);
732: return(0);
733: }
737: PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal time,Vec T,void *ctx)
738: {
740: PetscScalar *array;
741: PetscInt itime=(PetscInt)(time);
742: const TSType type;
743: PetscBool sundials;
744: MonitorCtx *user = (MonitorCtx*)ctx;
745: PetscViewer viewer = user->drawviewer;
746: PetscMPIInt rank;
747: PetscReal norm;
750: MPI_Comm_rank(((PetscObject)ts)->comm,&rank);
752: TSGetType(ts,&type);
753: PetscTypeCompare((PetscObject)ts,TSSUNDIALS,&sundials);
754: VecNorm(T,NORM_INFINITY,&norm);
755: if (sundials || itime%60 == 0){
756: VecGetArray(T,&array);
757: if (!rank){printf("step %4d, time %8.1f, %6.4f, %6.4f, %6.4f, %6.4f, %6.4f, %6.4f\n",step,time,(((array[0]-273)*9)/5 + 32),(((array[1]-273)*9)/5 + 32),array[2],array[3],array[4],array[5]);}
758: VecRestoreArray(T,&array);
760: if (user->drawcontours){
761: VecView(T,viewer);
762: }
763: }
764: return(0);
765: }