Actual source code: ex5.c

petsc-3.3-p6 2013-02-11
  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
122: extern PetscScalar emission(PetscScalar);                           //sets emission/absorption constant depending on water vapor content
123: extern PetscScalar calc_q(PetscScalar);                             //calculates specific humidity
124: extern PetscScalar mph2mpers(PetscScalar);                          //converts miles per hour to meters per second
125: 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.
126: extern PetscScalar fahr_to_cel(PetscScalar);                        //converts Fahrenheit to Celsius
127: extern PetscScalar cel_to_fahr(PetscScalar);                        //converts Celsius to Fahrenheit
128: extern PetscScalar calcmixingr(PetscScalar, PetscScalar);           //calculates mixing ratio
129: extern PetscScalar cloud(PetscScalar);                              //cloud radiative parameterization
130: extern PetscErrorCode FormInitialSolution(DM,Vec,void*);            //Specifies initial conditions for the system of equations (PETSc defined function)
131: extern PetscErrorCode RhsFunc(TS,PetscReal,Vec,Vec,void*);            //Specifies the user defined functions                       (PETSc defined function)
132: extern PetscErrorCode Monitor(TS,PetscInt,PetscReal,Vec,void*);            //Specifies output and visualization tools                       (PETSc defined function)
133: extern void readinput(struct in *put);                                    //reads input from text file
134: extern PetscErrorCode calcfluxs(PetscScalar, PetscScalar, PetscScalar, PetscScalar, PetscScalar, PetscScalar*);//calculates upward IR from surface
135: extern PetscErrorCode calcfluxa(PetscScalar, PetscScalar, PetscScalar, PetscScalar*);                          //calculates downward IR from atmosphere
136: extern PetscErrorCode sensibleflux(PetscScalar, PetscScalar, PetscScalar, PetscScalar*);                       //calculates sensible heat flux
137: extern PetscErrorCode potential_temperature(PetscScalar, PetscScalar, PetscScalar, PetscScalar, PetscScalar*); //calculates potential temperature
138: extern PetscErrorCode latentflux(PetscScalar, PetscScalar, PetscScalar, PetscScalar, PetscScalar*);            //calculates latent heat flux
139: extern PetscErrorCode calc_gflux(PetscScalar, PetscScalar, PetscScalar*);                                      //calculates flux between top soil layer and underlying earth

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;
163:   TS          ts;
164:   SNES        snes;
165:   DM          da;
166:   Vec         T,rhs;            /* solution vector */
167:   Mat         J;                /* Jacobian matrix */
168:   PetscReal   ftime,dt;
169:   PetscInt    steps,dof = 5;

171:   PetscInitialize(&argc,&argv,(char *)0,help);
172:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
173:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

175:   /* Inputs */
176:   readinput(&put);

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

184:   if (!rank) PetscPrintf(PETSC_COMM_SELF,"Initial Temperature = %g\n",sfctemp); //input surface temperature

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

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

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

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

208:     if (!rank){printf("Initial RH = %.1f percent\n\n",rh);}   //prints initial relative humidity

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

212:   /* Configure PETSc TS solver */
213:   /*------------------------------------------*/

215:   /*Create grid*/
216:   DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_PERIODIC,DMDA_BOUNDARY_PERIODIC,DMDA_STENCIL_STAR,-20,-20,
217:                       PETSC_DECIDE,PETSC_DECIDE,dof,1,PETSC_NULL,PETSC_NULL,&da);

219:   /*Define output window for each variable of interest*/
220:   DMDASetFieldName(da,0,"Ts");
221:   DMDASetFieldName(da,1,"Ta");
222:   DMDASetFieldName(da,2,"u");
223:   DMDASetFieldName(da,3,"v");
224:   DMDASetFieldName(da,4,"p");

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

242:   /*set values for MonitorCtx*/
243:   usermonitor.drawcontours = PETSC_FALSE;
244:   PetscOptionsHasName(PETSC_NULL,"-drawcontours",&usermonitor.drawcontours);
245:   if (usermonitor.drawcontours){
246:     PetscReal bounds[] = {1000.0,-1000.,  -1000.,-1000.,  1000.,-1000.,  1000.,-1000.,  1000,-1000, 100700,100800};
247:     PetscViewerDrawOpen(PETSC_COMM_WORLD,0,0,0,0,300,300,&usermonitor.drawviewer);
248:     PetscViewerDrawSetBounds(usermonitor.drawviewer,dof,bounds);
249:   }

251:   /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
252:      Extract global vectors from DA; 
253:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
254:   DMCreateGlobalVector(da,&T);
255:   VecDuplicate(T,&rhs); //r: vector to put the computed right hand side

257:   TSCreate(PETSC_COMM_WORLD,&ts);
258:   TSSetProblemType(ts,TS_NONLINEAR);
259:   TSSetType(ts,TSBEULER);
260:   TSSetRHSFunction(ts,rhs,RhsFunc,&user);

262:   /* Set Jacobian evaluation routine - use coloring to compute finite difference Jacobian efficiently */
263:   PetscBool      use_coloring=PETSC_TRUE;
264:   MatFDColoring  matfdcoloring=0;
265:   DMCreateMatrix(da,MATAIJ,&J);
266:   TSGetSNES(ts,&snes);
267:   if (use_coloring){
268:     ISColoring     iscoloring;
269:     DMCreateColoring(da,IS_COLORING_GLOBAL,MATAIJ,&iscoloring);
270:     MatFDColoringCreate(J,iscoloring,&matfdcoloring);
271:     MatFDColoringSetFromOptions(matfdcoloring);
272:     ISColoringDestroy(&iscoloring);
273:     MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))SNESTSFormFunction,ts);
274:     SNESSetJacobian(snes,J,J,SNESDefaultComputeJacobianColor,matfdcoloring);
275:   } else {
276:     SNESSetJacobian(snes,J,J,SNESDefaultComputeJacobian,PETSC_NULL);
277:   }

279:   /*Define what to print for ts_monitor option*/
280:   TSMonitorSet(ts,Monitor,&usermonitor,PETSC_NULL);
281:   FormInitialSolution(da,T,&user);
282:   dt    = TIMESTEP; /* initial time step */
283:   ftime = TIMESTEP*time;
284:   if (!rank){printf("time %d, ftime %g hour, TIMESTEP %g\n",time,ftime/3600,dt);}
285:   TSSetInitialTimeStep(ts,0.0,dt);
286:   TSSetDuration(ts,time,ftime);
287:   TSSetSolution(ts,T);

289:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
290:      Set runtime options
291:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
292:   TSSetFromOptions(ts);

294:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
295:      Solve nonlinear system
296:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
297:   TSSolve(ts,T,&ftime);
298:   TSGetTimeStepNumber(ts,&steps);
299:   if (!rank){PetscPrintf(PETSC_COMM_WORLD,"Solution T after %g hours %d steps\n",ftime/3600,steps);}


302:   if (matfdcoloring){MatFDColoringDestroy(&matfdcoloring);}
303:   if (usermonitor.drawcontours){
304:     PetscViewerDestroy(&usermonitor.drawviewer);
305:   }
306:   MatDestroy(&J);
307:   VecDestroy(&T);
308:   VecDestroy(&rhs);
309:   TSDestroy(&ts);
310:   DMDestroy(&da);

312:   PetscFinalize();
313:   return 0;
314: }
315: /*****************************end main program********************************/
316: /*****************************************************************************/
317: /*****************************************************************************/
318: /*****************************************************************************/
321: PetscErrorCode calcfluxs(PetscScalar sfctemp, PetscScalar airtemp, PetscScalar emma, PetscScalar fract, PetscScalar cloudTemp, PetscScalar* flux)
322: {
324:   *flux = SIG*((EMMSFC*emma*pow(airtemp,4)) + (EMMSFC*fract*(1 - emma)*pow(cloudTemp,4)) - (EMMSFC*pow(sfctemp,4)));   //calculates flux using Stefan-Boltzmann relation
325:   return(0);
326: }

330: PetscErrorCode calcfluxa(PetscScalar sfctemp, PetscScalar airtemp, PetscScalar emma, PetscScalar* flux)   //this function is not currently called upon
331: {
332:     PetscScalar emm = 0.001;
334:     *flux = SIG*(- emm*(pow(airtemp,4)));     //calculates flux usinge Stefan-Boltzmann relation
335:     return(0);
336: }
339: PetscErrorCode sensibleflux(PetscScalar sfctemp, PetscScalar airtemp, PetscScalar wind, PetscScalar* sheat)
340: {
341:       PetscScalar density = 1; //air density
342:       PetscScalar Cp = 1005;   //heat capicity for dry air
343:       PetscScalar wndmix;      //temperature change from wind mixing: wind*Ch


347:       wndmix = 0.0025 + 0.0042*wind;                               //regression equation valid for neutral and stable BL
348:       *sheat = density*Cp*wndmix*(airtemp - sfctemp);               //calculates sensible heat flux

350:       return(0);
351: }

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

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

376:      *latentheat = density*wndmix*beta*lhcnst*(q - qs); //calculates latent heat flux
377:      return(0);
378: }

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

390:      mixratio = calcmixingr(sfctemp,pressure1);

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

396:      pavg = ((0.7*pressure1)+pressure2)/2;         //calculates simple average press
397:      *pottemp = temp*(pow((pressure1/pavg),kdry)); //calculates potential temperature

399:      return(0);
400: }
401: extern PetscScalar calcmixingr(PetscScalar dtemp, PetscScalar pressure1)
402: {
403:      PetscScalar e;        //vapor pressure
404:      PetscScalar mixratio; //mixing ratio

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

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

423: PetscErrorCode calc_gflux(PetscScalar sfctemp, PetscScalar deep_grnd_temp, PetscScalar* Gflux)
424: {
425:       PetscScalar k;                       //thermal conductivity parameter
426:       PetscScalar n = 0.38;                   //value of soil porosity
427:       PetscScalar dz = 1;                  //depth of layer between soil surface and deep soil layer
428:       PetscScalar unit_soil_weight = 2700; //unit soil weight in kg/m^3


432:       k = ((0.135*(1-n)*unit_soil_weight) + 64.7)/(unit_soil_weight - (0.947*(1-n)*unit_soil_weight));  //dry soil conductivity
433:       *Gflux = (k*(deep_grnd_temp - sfctemp)/dz);   //calculates flux from deep ground layer
434:       return(0);
435: }
438: extern PetscScalar emission(PetscScalar pwat)
439: {
440:     PetscScalar emma;

442:     emma = 0.725 + 0.17*log10(pwat);

444:     return emma;
445: }
446: extern PetscScalar cloud(PetscScalar fract)
447: {
448:     PetscScalar emma = 0;

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

493:         ifp = fopen("ex5_control.txt", "r");

495:         for (i=0;i<110;i++)
496:                 fscanf(ifp, "%c", &x);
497:         fscanf(ifp, "%lf", &put->Ts);

499:         for (i=0;i<43;i++)
500:                 fscanf(ifp, "%c", &x);
501:         fscanf(ifp, "%lf", &put->Td);

503:         for (i=0;i<43;i++)
504:                 fscanf(ifp, "%c", &x);
505:         fscanf(ifp, "%lf", &put->Ta);

507:         for (i=0;i<43;i++)
508:                 fscanf(ifp, "%c", &x);
509:         fscanf(ifp, "%lf", &put->Tc);

511:         for (i=0;i<43;i++)
512:                 fscanf(ifp, "%c", &x);
513:         fscanf(ifp, "%lf", &put->fr);

515:         for (i=0;i<43;i++)
516:                 fscanf(ifp, "%c", &x);
517:         fscanf(ifp, "%lf", &put->wnd);

519:         for (i=0;i<43;i++)
520:                 fscanf(ifp, "%c", &x);
521:         fscanf(ifp, "%lf", &put->pwt);

523:         for (i=0;i<43;i++)
524:                 fscanf(ifp, "%c", &x);
525:         fscanf(ifp, "%lf", &put->wndDir);

527:         for (i=0;i<43;i++)
528:                 fscanf(ifp, "%c", &x);
529:         fscanf(ifp, "%lf", &put->time);

531:         for (i=0;i<63;i++)
532:                 fscanf(ifp, "%c", &x);
533:         fscanf(ifp, "%lf", &put->init);
534: }

536: /* ------------------------------------------------------------------- */
539: PetscErrorCode FormInitialSolution(DM da,Vec Xglobal,void *ctx)
540: {
542:   AppCtx         *user = (AppCtx*)ctx;       /* user-defined application context */
543:   PetscInt       i,j,xs,ys,xm,ym,Mx,My;
544:   Field          **X;
545:   PetscScalar    deltT;
546:   PetscReal      hx,hy;
547:   FILE           *ifp;
548:   FILE           *ofp;

551:   ofp = fopen("swing", "w");
552:   ifp = fopen("grid.in", "r");
553:   deltT = 0.8;

555:   DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
556:                      PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

558:   hx = 1/(PetscReal)(Mx-1);
559:   hy = 1/(PetscReal)(My-1);

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

564:   /* Get local grid boundaries */
565:   DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);

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

569:   if (user->init == 1) {
570:     for (j=ys; j<ys+ym; j++) {
571:       for (i=xs; i<xs+xm; i++) {
572:           X[j][i].Ts = user->Ts - i*0.0001;
573:           X[j][i].Ta = X[j][i].Ts - 5;
574:           X[j][i].u  = 0;
575:           X[j][i].v  = 0;
576:           X[j][i].p  = 1.25;
577:           if ((j == 5 || j == 6) && (i == 4 || i == 5)) X[j][i].p  += 0.00001;
578:           if ((j == 5 || j == 6) && (i == 12 || i == 13)) X[j][i].p  += 0.00001;
579:        }
580:     }
581:   }

583:   else {
584:     for (j=ys; j<ys+ym; j++) {
585:       for (i=xs; i<xs+xm; i++) {
586:           X[j][i].Ts = user->Ts;
587:           X[j][i].Ta = X[j][i].Ts - 5;
588:           X[j][i].u  = 0;
589:           X[j][i].v  = 0;
590:           X[j][i].p  = 1.25;
591:       }
592:     }
593:   }

595:   /* Restore vectors */
596:   DMDAVecRestoreArray(da,Xglobal,&X);
597:   return(0);
598: }

602: /*
603:    RhsFunc - Evaluates nonlinear function F(u).

605:    Input Parameters:
606: .  ts - the TS context
607: .  t - current time
608: .  Xglobal - input vector
609: .  F - output vector
610: .  ptr - optional user-defined context, as set by SNESSetFunction()

612:    Output Parameter:
613: .  F - rhs function vector
614:  */
615: PetscErrorCode RhsFunc(TS ts,PetscReal t,Vec Xglobal,Vec F,void *ctx)
616: {
617:   AppCtx         *user = (AppCtx*)ctx;       /* user-defined application context */
618:   DM             da = user->da;
620:   PetscInt       i,j,Mx,My,xs,ys,xm,ym;
621:   PetscReal      dhx,dhy;
622:   Vec            localT;
623:   Field          **X,**Frhs;                      //structures that contain variables of interest and left hand side of governing equations respectively
624:   PetscScalar    csoil     = user->csoil;     //heat constant for layer
625:   PetscScalar    dzlay     = user->dzlay;     //thickness of top soil layer
626:   PetscScalar    emma      = user->emma;      //emission parameter
627:   PetscScalar    wind      = user->wind;      //wind speed
628:   PetscScalar    dewtemp   = user->dewtemp;   //dew point temperature (moisture in air)
629:   PetscScalar    pressure1 = user->pressure1; //sea level pressure
630:   PetscScalar    airtemp   = user->airtemp;   //temperature of air near boundary layer inversion
631:   PetscScalar    fract     = user->fract;     //fraction of the sky covered by clouds
632:   PetscScalar    Tc           = user->Tc;              //temperature at base of lowest cloud layer
633:   PetscScalar    lat           = user->lat;              //latitude
634:   PetscScalar    Cp           = 1005.7;              //specific heat of air at constant pressure
635:   PetscScalar    Rd           = 287.058;              //gas constant for dry air
636:   PetscScalar    diffconst = 1000;              //diffusion coefficient
637:   PetscScalar    f           = 2*0.0000727*sin(lat);      //coriolis force
638:   PetscScalar    deep_grnd_temp = user->deep_grnd_temp; //temp in lowest ground layer
639:   PetscScalar    Ts,u,v,p,P;
640:   PetscScalar    u_abs,u_plus,u_minus,v_abs,v_plus,v_minus;

642:   PetscScalar         sfctemp1,fsfc1,Ra;
643:   PetscScalar         sheat;           //sensible heat flux
644:   PetscScalar         latentheat;      //latent heat flux
645:   PetscScalar         groundflux;      //flux from conduction of deep ground layer in contact with top soil
646:   PetscInt       xend,yend;

649:   DMGetLocalVector(da,&localT);
650:   DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
651:                      PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

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


657:   /*
658:      Scatter ghost points to local vector,using the 2-step process
659:         DAGlobalToLocalBegin(),DAGlobalToLocalEnd().
660:      By placing code between these two statements, computations can be
661:      done while messages are in transition.
662:   */
663:   DMGlobalToLocalBegin(da,Xglobal,INSERT_VALUES,localT);
664:   DMGlobalToLocalEnd(da,Xglobal,INSERT_VALUES,localT);

666:   /* Get pointers to vector data */
667:   DMDAVecGetArray(da,localT,&X);
668:   DMDAVecGetArray(da,F,&Frhs);

670:   /* Get local grid boundaries */
671:   DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);

673:   /* Compute function over the locally owned part of the grid */
674:   /* the interior points */
675:   xend=xs+xm; yend=ys+ym;
676:   for (j=ys; j<yend; j++) {
677:     for (i=xs; i<xend; i++) {
678:       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;

680:       sfctemp1 = (double)Ts;
681:       sfctemp1 = (double)X[j][i].Ts;
682:       calcfluxs(sfctemp1,airtemp,emma,fract,Tc,&fsfc1);       //calculates surface net radiative flux
683:       sensibleflux(sfctemp1,airtemp,wind,&sheat);             //calculate sensible heat flux
684:       latentflux(sfctemp1,dewtemp,wind,pressure1,&latentheat);//calculates latent heat flux
685:       calc_gflux(sfctemp1,deep_grnd_temp,&groundflux);        //calculates flux from earth below surface soil layer by conduction
686:       calcfluxa(sfctemp1,airtemp,emma,&Ra);                                         //Calculates the change in downward radiative flux
687:       fsfc1 = fsfc1 + latentheat + sheat + groundflux;                                         //adds radiative, sensible heat, latent heat, and ground heat flux yielding net flux

689:       /* convective coefficients for upwinding */
690:       u_abs = PetscAbsScalar(u);
691:       u_plus  = .5*(u + u_abs); // u if u>0; 0 if u<0
692:       u_minus = .5*(u - u_abs); // u if u <0; 0 if u>0

694:       v_abs = PetscAbsScalar(v);
695:       v_plus  = .5*(v + v_abs); // v if v>0; 0 if v<0
696:       v_minus = .5*(v - v_abs); // v if v <0; 0 if v>0

698:       /* Solve governing equations */
699:       P = p*Rd*Ts;

701:       /* du/dt -> time change of east-west component of the wind */
702:       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)
703:                         - v_plus*(u - X[j-1][i].u)*dhy - v_minus*(X[j+1][i].u - u)*dhy        // - v(du/dy)
704:                               -(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)]
705: //                        -(1/p)*(X[j][i+1].P - X[j][i-1].P)*dhx
706:                         + f*v;

708:       /* dv/dt -> time change of north-south component of the wind */
709:       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)
710:                         - v_plus*(v - X[j-1][i].v)*dhy - v_minus*(X[j+1][i].v - v)*dhy        // - v(dv/dy)
711:                               -(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)]
712: //                        -(1/p)*(X[j+1][i].P - X[j-1][i].P)*dhy
713:                         -f*u;

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

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

726:       Frhs[j][i].Ta = Ra/Cp;        //dTa/dt time change of air temperature
727:     }
728:   }

730:   /* Restore vectors */
731:   DMDAVecRestoreArray(da,localT,&X);
732:   DMDAVecRestoreArray(da,F,&Frhs);
733:   DMRestoreLocalVector(da,&localT);
734:   return(0);
735: }

739: PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal time,Vec T,void *ctx)
740: {
742:   PetscScalar    *array;
743:   PetscInt       itime=(PetscInt)(time);
744:   const TSType   type;
745:   PetscBool      sundials;
746:   MonitorCtx     *user = (MonitorCtx*)ctx;
747:   PetscViewer    viewer = user->drawviewer;
748:   PetscMPIInt    rank;
749:   PetscReal      norm;

752:   MPI_Comm_rank(((PetscObject)ts)->comm,&rank);

754:   TSGetType(ts,&type);
755:   PetscObjectTypeCompare((PetscObject)ts,TSSUNDIALS,&sundials);
756:   VecNorm(T,NORM_INFINITY,&norm);
757:   if (sundials || itime%60 == 0){
758:     VecGetArray(T,&array);
759:     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]);}
760:     VecRestoreArray(T,&array);

762:     if (user->drawcontours){
763:       VecView(T,viewer);
764:     }
765:   }
766:   return(0);
767: }