Actual source code: ex5.c

  1: static char help[] = "Nonlinear, time-dependent. Developed from radiative_surface_balance.c \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_view
  9:     ./ex5 -da_grid_x 20 -da_grid_y 20 -log_view
 10:     ./ex5 -da_grid_x 20 -da_grid_y 20 -ts_type rosw -ts_atol 1.e-6 -ts_rtol 1.e-6
 11:     ./ex5 -drawcontours -draw_pause 0.1 -draw_fields 0,1,2,3,4
 12: */

 14: /*
 15:    -----------------------------------------------------------------------

 17:    Governing equations:

 19:         R      = s*(Ea*Ta^4 - Es*Ts^4)
 20:         SH     = p*Cp*Ch*wind*(Ta - Ts)
 21:         LH     = p*L*Ch*wind*B(q(Ta) - q(Ts))
 22:         G      = k*(Tgnd - Ts)/dz

 24:         Fnet   = R + SH + LH + G

 26:         du/dt  = -u*(du/dx) - v*(du/dy) - 2*omeg*sin(lat)*v - (1/p)*(dP/dx)
 27:         dv/dt  = -u*(dv/dx) - v*(dv/dy) + 2*omeg*sin(lat)*u - (1/p)*(dP/dy)
 28:         dTs/dt = Fnet/(Cp*dz) - Div([u*Ts, v*Ts]) + D*Lap(Ts)
 29:                = Fnet/(Cs*dz) - u*(dTs/dx) - v*(dTs/dy) + D*(Ts_xx + Ts_yy)
 30:         dp/dt  = -Div([u*p,v*p])
 31:                = - u*dp/dx - v*dp/dy
 32:         dTa/dt = Fnet/Cp

 34:    Equation of State:

 36:         P = p*R*Ts

 38:    -----------------------------------------------------------------------

 40:    Program considers the evolution of a two dimensional atmosphere from
 41:    sunset to sunrise. There are two components:
 42:                 1. Surface energy balance model to compute diabatic dT (Fnet)
 43:                 2. Dynamical model using simplified primitive equations

 45:    Program is to be initiated at sunset and run to sunrise.

 47:    Inputs are:
 48:                 Surface temperature
 49:                 Dew point temperature
 50:                 Air temperature
 51:                 Temperature at cloud base (if clouds are present)
 52:                 Fraction of sky covered by clouds
 53:                 Wind speed
 54:                 Precipitable water in centimeters
 55:                 Wind direction

 57:    Inputs are are read in from the text file ex5_control.txt. To change an
 58:    input value use ex5_control.txt.

 60:    Solvers:
 61:             Backward Euler = default solver
 62:             Sundials = fastest and most accurate, requires Sundials libraries

 64:    This model is under development and should be used only as an example
 65:    and not as a predictive weather model.
 66: */

 68: #include <petscts.h>
 69: #include <petscdm.h>
 70: #include <petscdmda.h>

 72: /* stefan-boltzmann constant */
 73: #define SIG 0.000000056703
 74: /* absorption-emission constant for surface */
 75: #define EMMSFC   1
 76: /* amount of time (seconds) that passes before new flux is calculated */
 77: #define TIMESTEP 1

 79: /* variables of interest to be solved at each grid point */
 80: typedef struct {
 81:   PetscScalar Ts,Ta; /* surface and air temperature */
 82:   PetscScalar u,v;   /* wind speed */
 83:   PetscScalar p;     /* density */
 84: } Field;

 86: /* User defined variables. Used in solving for variables of interest */
 87: typedef struct {
 88:   DM          da;        /* grid */
 89:   PetscScalar csoil;     /* heat constant for layer */
 90:   PetscScalar dzlay;     /* thickness of top soil layer */
 91:   PetscScalar emma;      /* emission parameter */
 92:   PetscScalar wind;      /* wind speed */
 93:   PetscScalar dewtemp;   /* dew point temperature (moisture in air) */
 94:   PetscScalar pressure1; /* sea level pressure */
 95:   PetscScalar airtemp;   /* temperature of air near boundary layer inversion */
 96:   PetscScalar Ts;        /* temperature at the surface */
 97:   PetscScalar fract;     /* fraction of sky covered by clouds */
 98:   PetscScalar Tc;        /* temperature at base of lowest cloud layer */
 99:   PetscScalar lat;       /* Latitude in degrees */
100:   PetscScalar init;      /* initialization scenario */
101:   PetscScalar deep_grnd_temp; /* temperature of ground under top soil surface layer */
102: } AppCtx;

104: /* Struct for visualization */
105: typedef struct {
106:   PetscBool   drawcontours;   /* flag - 1 indicates drawing contours */
107:   PetscViewer drawviewer;
108:   PetscInt    interval;
109: } MonitorCtx;

111: /* Inputs read in from text file */
112: struct in {
113:   PetscScalar Ts;     /* surface temperature  */
114:   PetscScalar Td;     /* dewpoint temperature */
115:   PetscScalar Tc;     /* temperature of cloud base */
116:   PetscScalar fr;     /* fraction of sky covered by clouds */
117:   PetscScalar wnd;    /* wind speed */
118:   PetscScalar Ta;     /* air temperature */
119:   PetscScalar pwt;    /* precipitable water */
120:   PetscScalar wndDir; /* wind direction */
121:   PetscScalar lat;    /* latitude */
122:   PetscReal   time;   /* time in hours */
123:   PetscScalar init;
124: };

126: /* functions */
127: extern PetscScalar emission(PetscScalar);                           /* sets emission/absorption constant depending on water vapor content */
128: extern PetscScalar calc_q(PetscScalar);                             /* calculates specific humidity */
129: extern PetscScalar mph2mpers(PetscScalar);                          /* converts miles per hour to meters per second */
130: extern PetscScalar Lconst(PetscScalar);                             /* calculates latent heat constant taken from Satellite estimates of wind speed and latent heat flux over the global oceans., Bentamy et al. */
131: extern PetscScalar fahr_to_cel(PetscScalar);                        /* converts Fahrenheit to Celsius */
132: extern PetscScalar cel_to_fahr(PetscScalar);                        /* converts Celsius to Fahrenheit */
133: extern PetscScalar calcmixingr(PetscScalar, PetscScalar);           /* calculates mixing ratio */
134: extern PetscScalar cloud(PetscScalar);                              /* cloud radiative parameterization */
135: extern PetscErrorCode FormInitialSolution(DM,Vec,void*);            /* Specifies initial conditions for the system of equations (PETSc defined function) */
136: extern PetscErrorCode RhsFunc(TS,PetscReal,Vec,Vec,void*);          /* Specifies the user defined functions                     (PETSc defined function) */
137: extern PetscErrorCode Monitor(TS,PetscInt,PetscReal,Vec,void*);     /* Specifies output and visualization tools                 (PETSc defined function) */
138: extern PetscErrorCode readinput(struct in *put);                              /* reads input from text file */
139: extern PetscErrorCode calcfluxs(PetscScalar, PetscScalar, PetscScalar, PetscScalar, PetscScalar, PetscScalar*); /* calculates upward IR from surface */
140: extern PetscErrorCode calcfluxa(PetscScalar, PetscScalar, PetscScalar, PetscScalar*);                           /* calculates downward IR from atmosphere */
141: extern PetscErrorCode sensibleflux(PetscScalar, PetscScalar, PetscScalar, PetscScalar*);                        /* calculates sensible heat flux */
142: extern PetscErrorCode potential_temperature(PetscScalar, PetscScalar, PetscScalar, PetscScalar, PetscScalar*);  /* calculates potential temperature */
143: extern PetscErrorCode latentflux(PetscScalar, PetscScalar, PetscScalar, PetscScalar, PetscScalar*);             /* calculates latent heat flux */
144: extern PetscErrorCode calc_gflux(PetscScalar, PetscScalar, PetscScalar*);                                       /* calculates flux between top soil layer and underlying earth */

146: int main(int argc,char **argv)
147: {
149:   PetscInt       time;           /* amount of loops */
150:   struct in      put;
151:   PetscScalar    rh;             /* relative humidity */
152:   PetscScalar    x;              /* memory varialbe for relative humidity calculation */
153:   PetscScalar    deep_grnd_temp; /* temperature of ground under top soil surface layer */
154:   PetscScalar    emma;           /* absorption-emission constant for air */
155:   PetscScalar    pressure1 = 101300; /* surface pressure */
156:   PetscScalar    mixratio;       /* mixing ratio */
157:   PetscScalar    airtemp;        /* temperature of air near boundary layer inversion */
158:   PetscScalar    dewtemp;        /* dew point temperature */
159:   PetscScalar    sfctemp;        /* temperature at surface */
160:   PetscScalar    pwat;           /* total column precipitable water */
161:   PetscScalar    cloudTemp;      /* temperature at base of cloud */
162:   AppCtx         user;           /*  user-defined work context */
163:   MonitorCtx     usermonitor;    /* user-defined monitor context */
164:   TS             ts;
165:   SNES           snes;
166:   DM             da;
167:   Vec            T,rhs;          /* solution vector */
168:   Mat            J;              /* Jacobian matrix */
169:   PetscReal      ftime,dt;
170:   PetscInt       steps,dof = 5;
171:   PetscBool      use_coloring  = PETSC_TRUE;
172:   MatFDColoring  matfdcoloring = 0;
173:   PetscBool      monitor_off = PETSC_FALSE;

175:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;

177:   /* Inputs */
178:   readinput(&put);

180:   sfctemp   = put.Ts;
181:   dewtemp   = put.Td;
182:   cloudTemp = put.Tc;
183:   airtemp   = put.Ta;
184:   pwat      = put.pwt;

186:   PetscPrintf(PETSC_COMM_WORLD,"Initial Temperature = %g\n",(double)sfctemp); /* input surface temperature */

188:   deep_grnd_temp = sfctemp - 10;   /* set underlying ground layer temperature */
189:   emma           = emission(pwat); /* accounts for radiative effects of water vapor */

191:   /* Converts from Fahrenheit to Celsuis */
192:   sfctemp        = fahr_to_cel(sfctemp);
193:   airtemp        = fahr_to_cel(airtemp);
194:   dewtemp        = fahr_to_cel(dewtemp);
195:   cloudTemp      = fahr_to_cel(cloudTemp);
196:   deep_grnd_temp = fahr_to_cel(deep_grnd_temp);

198:   /* Converts from Celsius to Kelvin */
199:   sfctemp        += 273;
200:   airtemp        += 273;
201:   dewtemp        += 273;
202:   cloudTemp      += 273;
203:   deep_grnd_temp += 273;

205:   /* Calculates initial relative humidity */
206:   x        = calcmixingr(dewtemp,pressure1);
207:   mixratio = calcmixingr(sfctemp,pressure1);
208:   rh       = (x/mixratio)*100;

210:   PetscPrintf(PETSC_COMM_WORLD,"Initial RH = %.1f percent\n\n",(double)rh);   /* prints initial relative humidity */

212:   time = 3600*put.time;                         /* sets amount of timesteps to run model */

214:   /* Configure PETSc TS solver */
215:   /*------------------------------------------*/

217:   /* Create grid */
218:   DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,DM_BOUNDARY_PERIODIC,DMDA_STENCIL_STAR,20,20,PETSC_DECIDE,PETSC_DECIDE,dof,1,NULL,NULL,&da);
219:   DMSetFromOptions(da);
220:   DMSetUp(da);
221:   DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);

223:   /* Define output window for each variable of interest */
224:   DMDASetFieldName(da,0,"Ts");
225:   DMDASetFieldName(da,1,"Ta");
226:   DMDASetFieldName(da,2,"u");
227:   DMDASetFieldName(da,3,"v");
228:   DMDASetFieldName(da,4,"p");

230:   /* set values for appctx */
231:   user.da             = da;
232:   user.Ts             = sfctemp;
233:   user.fract          = put.fr;          /* fraction of sky covered by clouds */
234:   user.dewtemp        = dewtemp;         /* dew point temperature (mositure in air) */
235:   user.csoil          = 2000000;         /* heat constant for layer */
236:   user.dzlay          = 0.08;            /* thickness of top soil layer */
237:   user.emma           = emma;            /* emission parameter */
238:   user.wind           = put.wnd;         /* wind spped */
239:   user.pressure1      = pressure1;       /* sea level pressure */
240:   user.airtemp        = airtemp;         /* temperature of air near boundar layer inversion */
241:   user.Tc             = cloudTemp;       /* temperature at base of lowest cloud layer */
242:   user.init           = put.init;        /* user chosen initiation scenario */
243:   user.lat            = 70*0.0174532;    /* converts latitude degrees to latitude in radians */
244:   user.deep_grnd_temp = deep_grnd_temp;  /* temp in lowest ground layer */

246:   /* set values for MonitorCtx */
247:   usermonitor.drawcontours = PETSC_FALSE;
248:   PetscOptionsHasName(NULL,NULL,"-drawcontours",&usermonitor.drawcontours);
249:   if (usermonitor.drawcontours) {
250:     PetscReal bounds[] = {1000.0,-1000.,  -1000.,-1000.,  1000.,-1000.,  1000.,-1000.,  1000,-1000, 100700,100800};
251:     PetscViewerDrawOpen(PETSC_COMM_WORLD,0,0,0,0,300,300,&usermonitor.drawviewer);
252:     PetscViewerDrawSetBounds(usermonitor.drawviewer,dof,bounds);
253:   }
254:   usermonitor.interval = 1;
255:   PetscOptionsGetInt(NULL,NULL,"-monitor_interval",&usermonitor.interval,NULL);

257:   /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
258:      Extract global vectors from DA;
259:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
260:   DMCreateGlobalVector(da,&T);
261:   VecDuplicate(T,&rhs); /* r: vector to put the computed right hand side */

263:   TSCreate(PETSC_COMM_WORLD,&ts);
264:   TSSetProblemType(ts,TS_NONLINEAR);
265:   TSSetType(ts,TSBEULER);
266:   TSSetRHSFunction(ts,rhs,RhsFunc,&user);

268:   /* Set Jacobian evaluation routine - use coloring to compute finite difference Jacobian efficiently */
269:   DMSetMatType(da,MATAIJ);
270:   DMCreateMatrix(da,&J);
271:   TSGetSNES(ts,&snes);
272:   if (use_coloring) {
273:     ISColoring iscoloring;
274:     DMCreateColoring(da,IS_COLORING_GLOBAL,&iscoloring);
275:     MatFDColoringCreate(J,iscoloring,&matfdcoloring);
276:     MatFDColoringSetFromOptions(matfdcoloring);
277:     MatFDColoringSetUp(J,iscoloring,matfdcoloring);
278:     ISColoringDestroy(&iscoloring);
279:     MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))SNESTSFormFunction,ts);
280:     SNESSetJacobian(snes,J,J,SNESComputeJacobianDefaultColor,matfdcoloring);
281:   } else {
282:     SNESSetJacobian(snes,J,J,SNESComputeJacobianDefault,NULL);
283:   }

285:   /* Define what to print for ts_monitor option */
286:   PetscOptionsHasName(NULL,NULL,"-monitor_off",&monitor_off);
287:   if (!monitor_off) {
288:     TSMonitorSet(ts,Monitor,&usermonitor,NULL);
289:   }
290:   FormInitialSolution(da,T,&user);
291:   dt    = TIMESTEP; /* initial time step */
292:   ftime = TIMESTEP*time;
293:   PetscPrintf(PETSC_COMM_WORLD,"time %D, ftime %g hour, TIMESTEP %g\n",time,(double)(ftime/3600),(double)dt);

295:   TSSetTimeStep(ts,dt);
296:   TSSetMaxSteps(ts,time);
297:   TSSetMaxTime(ts,ftime);
298:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
299:   TSSetSolution(ts,T);
300:   TSSetDM(ts,da);

302:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
303:      Set runtime options
304:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
305:   TSSetFromOptions(ts);

307:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
308:      Solve nonlinear system
309:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
310:   TSSolve(ts,T);
311:   TSGetSolveTime(ts,&ftime);
312:   TSGetStepNumber(ts,&steps);
313:   PetscPrintf(PETSC_COMM_WORLD,"Solution T after %g hours %D steps\n",(double)(ftime/3600),steps);

315:   if (matfdcoloring) {MatFDColoringDestroy(&matfdcoloring);}
316:   if (usermonitor.drawcontours) {
317:     PetscViewerDestroy(&usermonitor.drawviewer);
318:   }
319:   MatDestroy(&J);
320:   VecDestroy(&T);
321:   VecDestroy(&rhs);
322:   TSDestroy(&ts);
323:   DMDestroy(&da);

325:   PetscFinalize();
326:   return ierr;
327: }
328: /*****************************end main program********************************/
329: /*****************************************************************************/
330: /*****************************************************************************/
331: /*****************************************************************************/
332: PetscErrorCode calcfluxs(PetscScalar sfctemp, PetscScalar airtemp, PetscScalar emma, PetscScalar fract, PetscScalar cloudTemp, PetscScalar *flux)
333: {
335:   *flux = SIG*((EMMSFC*emma*PetscPowScalarInt(airtemp,4)) + (EMMSFC*fract*(1 - emma)*PetscPowScalarInt(cloudTemp,4)) - (EMMSFC*PetscPowScalarInt(sfctemp,4)));   /* calculates flux using Stefan-Boltzmann relation */
336:   return(0);
337: }

339: PetscErrorCode calcfluxa(PetscScalar sfctemp, PetscScalar airtemp, PetscScalar emma, PetscScalar *flux)   /* this function is not currently called upon */
340: {
341:   PetscScalar emm = 0.001;

344:   *flux = SIG*(-emm*(PetscPowScalarInt(airtemp,4)));      /* calculates flux usinge Stefan-Boltzmann relation */
345:   return(0);
346: }
347: PetscErrorCode sensibleflux(PetscScalar sfctemp, PetscScalar airtemp, PetscScalar wind, PetscScalar *sheat)
348: {
349:   PetscScalar density = 1;    /* air density */
350:   PetscScalar Cp      = 1005; /* heat capicity for dry air */
351:   PetscScalar wndmix;         /* temperature change from wind mixing: wind*Ch */

354:   wndmix = 0.0025 + 0.0042*wind;                               /* regression equation valid for neutral and stable BL */
355:   *sheat = density*Cp*wndmix*(airtemp - sfctemp);              /* calculates sensible heat flux */
356:   return(0);
357: }

359: PetscErrorCode latentflux(PetscScalar sfctemp, PetscScalar dewtemp, PetscScalar wind, PetscScalar pressure1, PetscScalar *latentheat)
360: {
361:   PetscScalar density = 1;   /* density of dry air */
362:   PetscScalar q;             /* actual specific humitity */
363:   PetscScalar qs;            /* saturation specific humidity */
364:   PetscScalar wndmix;        /* temperature change from wind mixing: wind*Ch */
365:   PetscScalar beta = .4;     /* moisture availability */
366:   PetscScalar mr;            /* mixing ratio */
367:   PetscScalar lhcnst;        /* latent heat of vaporization constant = 2501000 J/kg at 0c */
368:                              /* latent heat of saturation const = 2834000 J/kg */
369:                              /* latent heat of fusion const = 333700 J/kg */

372:   wind   = mph2mpers(wind);                /* converts wind from mph to meters per second */
373:   wndmix = 0.0025 + 0.0042*wind;           /* regression equation valid for neutral BL */
374:   lhcnst = Lconst(sfctemp);                /* calculates latent heat of evaporation */
375:   mr     = calcmixingr(sfctemp,pressure1); /* calculates saturation mixing ratio */
376:   qs     = calc_q(mr);                     /* calculates saturation specific humidty */
377:   mr     = calcmixingr(dewtemp,pressure1); /* calculates mixing ratio */
378:   q      = calc_q(mr);                     /* calculates specific humidty */

380:   *latentheat = density*wndmix*beta*lhcnst*(q - qs); /* calculates latent heat flux */
381:   return(0);
382: }

384: PetscErrorCode potential_temperature(PetscScalar temp, PetscScalar pressure1, PetscScalar pressure2, PetscScalar sfctemp, PetscScalar *pottemp)
385: {
386:   PetscScalar kdry;     /* poisson constant for dry atmosphere */
387:   PetscScalar pavg;     /* average atmospheric pressure */
388:   /* PetscScalar mixratio; mixing ratio */
389:   /* PetscScalar kmoist;   poisson constant for moist atmosphere */

392:   /* mixratio = calcmixingr(sfctemp,pressure1); */

394: /* initialize poisson constant */
395:   kdry   = 0.2854;
396:   /* kmoist = 0.2854*(1 - 0.24*mixratio); */

398:   pavg     = ((0.7*pressure1)+pressure2)/2;     /* calculates simple average press */
399:   *pottemp = temp*(PetscPowScalar((pressure1/pavg),kdry)); /* calculates potential temperature */
400:   return(0);
401: }
402: extern PetscScalar calcmixingr(PetscScalar dtemp, PetscScalar pressure1)
403: {
404:   PetscScalar e;        /* vapor pressure */
405:   PetscScalar mixratio; /* mixing ratio */

407:   dtemp    = dtemp - 273;                                /* converts from Kelvin to Celsuis */
408:   e        = 6.11*(PetscPowScalar(10,((7.5*dtemp)/(237.7+dtemp)))); /* converts from dew point temp to vapor pressure */
409:   e        = e*100;                                      /* converts from hPa to Pa */
410:   mixratio = (0.622*e)/(pressure1 - e);                  /* computes mixing ratio */
411:   mixratio = mixratio*1;                                 /* convert to g/Kg */

413:   return mixratio;
414: }
415: extern PetscScalar calc_q(PetscScalar rv)
416: {
417:   PetscScalar specific_humidity;        /* define specific humidity variable */
418:   specific_humidity = rv/(1 + rv);      /* calculates specific humidity */
419:   return specific_humidity;
420: }

422: PetscErrorCode calc_gflux(PetscScalar sfctemp, PetscScalar deep_grnd_temp, PetscScalar* Gflux)
423: {
424:   PetscScalar k;                       /* thermal conductivity parameter */
425:   PetscScalar n                = 0.38; /* value of soil porosity */
426:   PetscScalar dz               = 1;    /* depth of layer between soil surface and deep soil layer */
427:   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: }
434: extern PetscScalar emission(PetscScalar pwat)
435: {
436:   PetscScalar emma;

438:   emma = 0.725 + 0.17*PetscLog10Real(PetscRealPart(pwat));

440:   return emma;
441: }
442: extern PetscScalar cloud(PetscScalar fract)
443: {
444:   PetscScalar emma = 0;

446:   /* modifies radiative balance depending on cloud cover */
447:   if (fract >= 0.9)                     emma = 1;
448:   else if (0.9 > fract && fract >= 0.8) emma = 0.9;
449:   else if (0.8 > fract && fract >= 0.7) emma = 0.85;
450:   else if (0.7 > fract && fract >= 0.6) emma = 0.75;
451:   else if (0.6 > fract && fract >= 0.5) emma = 0.65;
452:   else if (0.4 > fract && fract >= 0.3) emma = emma*1.086956;
453:   return emma;
454: }
455: extern PetscScalar Lconst(PetscScalar sfctemp)
456: {
457:   PetscScalar Lheat;
458:   sfctemp -=273;                               /* converts from kelvin to celsius */
459:   Lheat    = 4186.8*(597.31 - 0.5625*sfctemp); /* calculates latent heat constant */
460:   return Lheat;
461: }
462: extern PetscScalar mph2mpers(PetscScalar wind)
463: {
464:   wind = ((wind*1.6*1000)/3600);                 /* converts wind from mph to meters per second */
465:   return wind;
466: }
467: extern PetscScalar fahr_to_cel(PetscScalar temp)
468: {
469:   temp = (5*(temp-32))/9; /* converts from farhrenheit to celsuis */
470:   return temp;
471: }
472: extern PetscScalar cel_to_fahr(PetscScalar temp)
473: {
474:   temp = ((temp*9)/5) + 32; /* converts from celsuis to farhrenheit */
475:   return temp;
476: }
477: PetscErrorCode readinput(struct in *put)
478: {
479:   int    i;
480:   char   x;
481:   FILE   *ifp;
482:   double tmp;

485:   ifp = fopen("ex5_control.txt", "r");
486:   if (!ifp) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_OPEN,"Unable to open input file");
487:   for (i=0; i<110; i++) { if (fscanf(ifp, "%c", &x) != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Unable to read file");}
488:   if (fscanf(ifp, "%lf", &tmp) != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Unable to read file");
489:   put->Ts = tmp;

491:   for (i=0; i<43; i++) { if (fscanf(ifp, "%c", &x) != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Unable to read file");}
492:   if (fscanf(ifp, "%lf", &tmp) != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Unable to read file");
493:   put->Td = tmp;

495:   for (i=0; i<43; i++) { if (fscanf(ifp, "%c", &x) != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Unable to read file");}
496:   if (fscanf(ifp, "%lf", &tmp) != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Unable to read file");
497:   put->Ta = tmp;

499:   for (i=0; i<43; i++) { if (fscanf(ifp, "%c", &x) != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Unable to read file");}
500:   if (fscanf(ifp, "%lf", &tmp)!= 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Unable to read file");
501:   put->Tc = tmp;

503:   for (i=0; i<43; i++) { if (fscanf(ifp, "%c", &x) != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Unable to read file");}
504:   if (fscanf(ifp, "%lf", &tmp) != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Unable to read file");
505:   put->fr = tmp;

507:   for (i=0; i<43; i++) {if (fscanf(ifp, "%c", &x) != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Unable to read file");}
508:   if (fscanf(ifp, "%lf", &tmp) != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Unable to read file");
509:   put->wnd = tmp;

511:   for (i=0; i<43; i++) {if (fscanf(ifp, "%c", &x) != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Unable to read file");}
512:   if (fscanf(ifp, "%lf", &tmp) != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Unable to read file");
513:   put->pwt = tmp;

515:   for (i=0; i<43; i++) {if (fscanf(ifp, "%c", &x) != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Unable to read file");}
516:   if (fscanf(ifp, "%lf", &tmp) != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Unable to read file");
517:   put->wndDir = tmp;

519:   for (i=0; i<43; i++) {if (fscanf(ifp, "%c", &x) != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Unable to read file");}
520:   if (fscanf(ifp, "%lf", &tmp) != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Unable to read file");
521:   put->time = tmp;

523:   for (i=0; i<63; i++) {if (fscanf(ifp, "%c", &x) != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Unable to read file");}
524:   if (fscanf(ifp, "%lf", &tmp) != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Unable to read file");
525:   put->init = tmp;
526:   return(0);
527: }

529: /* ------------------------------------------------------------------- */
530: PetscErrorCode FormInitialSolution(DM da,Vec Xglobal,void *ctx)
531: {
533:   AppCtx         *user = (AppCtx*)ctx;       /* user-defined application context */
534:   PetscInt       i,j,xs,ys,xm,ym,Mx,My;
535:   Field          **X;

538:   DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
539:                      PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

541:   /* Get pointers to vector data */
542:   DMDAVecGetArray(da,Xglobal,&X);

544:   /* Get local grid boundaries */
545:   DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);

547:   /* Compute function over the locally owned part of the grid */

549:   if (user->init == 1) {
550:     for (j=ys; j<ys+ym; j++) {
551:       for (i=xs; i<xs+xm; i++) {
552:         X[j][i].Ts = user->Ts - i*0.0001;
553:         X[j][i].Ta = X[j][i].Ts - 5;
554:         X[j][i].u  = 0;
555:         X[j][i].v  = 0;
556:         X[j][i].p  = 1.25;
557:         if ((j == 5 || j == 6) && (i == 4 || i == 5))   X[j][i].p += 0.00001;
558:         if ((j == 5 || j == 6) && (i == 12 || i == 13)) X[j][i].p += 0.00001;
559:       }
560:     }
561:   } else {
562:     for (j=ys; j<ys+ym; j++) {
563:       for (i=xs; i<xs+xm; i++) {
564:         X[j][i].Ts = user->Ts;
565:         X[j][i].Ta = X[j][i].Ts - 5;
566:         X[j][i].u  = 0;
567:         X[j][i].v  = 0;
568:         X[j][i].p  = 1.25;
569:       }
570:     }
571:   }

573:   /* Restore vectors */
574:   DMDAVecRestoreArray(da,Xglobal,&X);
575:   return(0);
576: }

578: /*
579:    RhsFunc - Evaluates nonlinear function F(u).

581:    Input Parameters:
582: .  ts - the TS context
583: .  t - current time
584: .  Xglobal - input vector
585: .  F - output vector
586: .  ptr - optional user-defined context, as set by SNESSetFunction()

588:    Output Parameter:
589: .  F - rhs function vector
590:  */
591: PetscErrorCode RhsFunc(TS ts,PetscReal t,Vec Xglobal,Vec F,void *ctx)
592: {
593:   AppCtx         *user = (AppCtx*)ctx;       /* user-defined application context */
594:   DM             da    = user->da;
596:   PetscInt       i,j,Mx,My,xs,ys,xm,ym;
597:   PetscReal      dhx,dhy;
598:   Vec            localT;
599:   Field          **X,**Frhs;                            /* structures that contain variables of interest and left hand side of governing equations respectively */
600:   PetscScalar    csoil          = user->csoil;          /* heat constant for layer */
601:   PetscScalar    dzlay          = user->dzlay;          /* thickness of top soil layer */
602:   PetscScalar    emma           = user->emma;           /* emission parameter */
603:   PetscScalar    wind           = user->wind;           /* wind speed */
604:   PetscScalar    dewtemp        = user->dewtemp;        /* dew point temperature (moisture in air) */
605:   PetscScalar    pressure1      = user->pressure1;      /* sea level pressure */
606:   PetscScalar    airtemp        = user->airtemp;        /* temperature of air near boundary layer inversion */
607:   PetscScalar    fract          = user->fract;          /* fraction of the sky covered by clouds */
608:   PetscScalar    Tc             = user->Tc;             /* temperature at base of lowest cloud layer */
609:   PetscScalar    lat            = user->lat;            /* latitude */
610:   PetscScalar    Cp             = 1005.7;               /* specific heat of air at constant pressure */
611:   PetscScalar    Rd             = 287.058;              /* gas constant for dry air */
612:   PetscScalar    diffconst      = 1000;                 /* diffusion coefficient */
613:   PetscScalar    f              = 2*0.0000727*PetscSinScalar(lat); /* coriolis force */
614:   PetscScalar    deep_grnd_temp = user->deep_grnd_temp; /* temp in lowest ground layer */
615:   PetscScalar    Ts,u,v,p;
616:   PetscScalar    u_abs,u_plus,u_minus,v_abs,v_plus,v_minus;

618:   PetscScalar sfctemp1,fsfc1,Ra;
619:   PetscScalar sheat;                   /* sensible heat flux */
620:   PetscScalar latentheat;              /* latent heat flux */
621:   PetscScalar groundflux;              /* flux from conduction of deep ground layer in contact with top soil */
622:   PetscInt    xend,yend;

625:   DMGetLocalVector(da,&localT);
626:   DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

628:   dhx = (PetscReal)(Mx-1)/(5000*(Mx-1));  /* dhx = 1/dx; assume 2D space domain: [0.0, 1.e5] x [0.0, 1.e5] */
629:   dhy = (PetscReal)(My-1)/(5000*(Mx-1));  /* dhy = 1/dy; */

631:   /*
632:      Scatter ghost points to local vector,using the 2-step process
633:         DAGlobalToLocalBegin(),DAGlobalToLocalEnd().
634:      By placing code between these two statements, computations can be
635:      done while messages are in transition.
636:   */
637:   DMGlobalToLocalBegin(da,Xglobal,INSERT_VALUES,localT);
638:   DMGlobalToLocalEnd(da,Xglobal,INSERT_VALUES,localT);

640:   /* Get pointers to vector data */
641:   DMDAVecGetArrayRead(da,localT,&X);
642:   DMDAVecGetArray(da,F,&Frhs);

644:   /* Get local grid boundaries */
645:   DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);

647:   /* Compute function over the locally owned part of the grid */
648:   /* the interior points */
649:   xend=xs+xm; yend=ys+ym;
650:   for (j=ys; j<yend; j++) {
651:     for (i=xs; i<xend; i++) {
652:       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; */

654:       sfctemp1 = (double)Ts;
655:       calcfluxs(sfctemp1,airtemp,emma,fract,Tc,&fsfc1);        /* calculates surface net radiative flux */
656:       sensibleflux(sfctemp1,airtemp,wind,&sheat);              /* calculate sensible heat flux */
657:       latentflux(sfctemp1,dewtemp,wind,pressure1,&latentheat); /* calculates latent heat flux */
658:       calc_gflux(sfctemp1,deep_grnd_temp,&groundflux);         /* calculates flux from earth below surface soil layer by conduction */
659:       calcfluxa(sfctemp1,airtemp,emma,&Ra);                    /* Calculates the change in downward radiative flux */
660:       fsfc1    = fsfc1 + latentheat + sheat + groundflux;                               /* adds radiative, sensible heat, latent heat, and ground heat flux yielding net flux */

662:       /* convective coefficients for upwinding */
663:       u_abs   = PetscAbsScalar(u);
664:       u_plus  = .5*(u + u_abs); /* u if u>0; 0 if u<0 */
665:       u_minus = .5*(u - u_abs); /* u if u <0; 0 if u>0 */

667:       v_abs   = PetscAbsScalar(v);
668:       v_plus  = .5*(v + v_abs); /* v if v>0; 0 if v<0 */
669:       v_minus = .5*(v - v_abs); /* v if v <0; 0 if v>0 */

671:       /* Solve governing equations */
672:       /* P = p*Rd*Ts; */

674:       /* du/dt -> time change of east-west component of the wind */
675:       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) */
676:                      - v_plus*(u - X[j-1][i].u)*dhy - v_minus*(X[j+1][i].u - u)*dhy       /* - v(du/dy) */
677:                      -(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)] */
678: /*                     -(1/p)*(X[j][i+1].P - X[j][i-1].P)*dhx */
679:                      + f*v;

681:       /* dv/dt -> time change of north-south component of the wind */
682:       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) */
683:                      - v_plus*(v - X[j-1][i].v)*dhy - v_minus*(X[j+1][i].v - v)*dhy       /* - v(dv/dy) */
684:                      -(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)] */
685: /*                     -(1/p)*(X[j+1][i].P - X[j-1][i].P)*dhy */
686:                      -f*u;

688:       /* dT/dt -> time change of temperature */
689:       Frhs[j][i].Ts = (fsfc1/(csoil*dzlay))                                            /* Fnet/(Cp*dz)  diabatic change in T */
690:                       -u_plus*(Ts - X[j][i-1].Ts)*dhx - u_minus*(X[j][i+1].Ts - Ts)*dhx  /* - u*(dTs/dx)  advection x */
691:                       -v_plus*(Ts - X[j-1][i].Ts)*dhy - v_minus*(X[j+1][i].Ts - Ts)*dhy  /* - v*(dTs/dy)  advection y */
692:                       + diffconst*((X[j][i+1].Ts - 2*Ts + X[j][i-1].Ts)*dhx*dhx               /* + D(Ts_xx + Ts_yy)  diffusion */
693:                                    + (X[j+1][i].Ts - 2*Ts + X[j-1][i].Ts)*dhy*dhy);

695:       /* dp/dt -> time change of */
696:       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) */
697:                      -v_plus*(p - X[j-1][i].p)*dhy - v_minus*(X[j+1][i].p - p)*dhy;    /* - v*(dp/dy) */

699:       Frhs[j][i].Ta = Ra/Cp;  /* dTa/dt time change of air temperature */
700:     }
701:   }

703:   /* Restore vectors */
704:   DMDAVecRestoreArrayRead(da,localT,&X);
705:   DMDAVecRestoreArray(da,F,&Frhs);
706:   DMRestoreLocalVector(da,&localT);
707:   return(0);
708: }

710: PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal time,Vec T,void *ctx)
711: {
712:   PetscErrorCode    ierr;
713:   const PetscScalar *array;
714:   MonitorCtx        *user  = (MonitorCtx*)ctx;
715:   PetscViewer       viewer = user->drawviewer;
716:   PetscReal         norm;

719:   VecNorm(T,NORM_INFINITY,&norm);

721:   if (step%user->interval == 0) {
722:     VecGetArrayRead(T,&array);
723:     PetscPrintf(PETSC_COMM_WORLD,"step %D, time %8.1f,  %6.4f, %6.4f, %6.4f, %6.4f, %6.4f, %6.4f\n",step,(double)time,(double)(((array[0]-273)*9)/5 + 32),(double)(((array[1]-273)*9)/5 + 32),(double)array[2],(double)array[3],(double)array[4],(double)array[5]);
724:     VecRestoreArrayRead(T,&array);
725:   }

727:   if (user->drawcontours) {
728:     VecView(T,viewer);
729:   }
730:   return(0);
731: }

733: /*TEST

735:    build:
736:       requires: !complex !single

738:    test:
739:       args: -ts_max_steps 130 -monitor_interval 60
740:       output_file: output/ex5.out
741:       requires: !complex !single
742:       localrunfiles: ex5_control.txt

744:    test:
745:       suffix: 2
746:       nsize: 4
747:       args: -ts_max_steps 130 -monitor_interval 60
748:       output_file: output/ex5.out
749:       localrunfiles: ex5_control.txt
750:       requires: !complex !single

752: TEST*/