1 static char help[] = "Nonlinear, time-dependent. Developed from radiative_surface_balance.c \n"; 2 /* 3 Contributed by Steve Froehlich, Illinois Institute of Technology 4 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 */ 13 14 /* 15 ----------------------------------------------------------------------- 16 17 Governing equations: 18 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 23 24 Fnet = R + SH + LH + G 25 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 33 34 Equation of State: 35 36 P = p*R*Ts 37 38 ----------------------------------------------------------------------- 39 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 44 45 Program is to be initiated at sunset and run to sunrise. 46 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 56 57 Inputs are are read in from the text file ex5_control.txt. To change an 58 input value use ex5_control.txt. 59 60 Solvers: 61 Backward Euler = default solver 62 Sundials = fastest and most accurate, requires Sundials libraries 63 64 This model is under development and should be used only as an example 65 and not as a predictive weather model. 66 */ 67 68 #include <petscts.h> 69 #include <petscdm.h> 70 #include <petscdmda.h> 71 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 78 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; 85 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; 103 104 /* Struct for visualization */ 105 typedef struct { 106 PetscBool drawcontours; /* flag - 1 indicates drawing contours */ 107 PetscViewer drawviewer; 108 PetscInt interval; 109 } MonitorCtx; 110 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 }; 125 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 */ 145 146 int main(int argc,char **argv) 147 { 148 PetscErrorCode ierr; 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; 174 175 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 176 177 /* Inputs */ 178 ierr = readinput(&put);CHKERRQ(ierr); 179 180 sfctemp = put.Ts; 181 dewtemp = put.Td; 182 cloudTemp = put.Tc; 183 airtemp = put.Ta; 184 pwat = put.pwt; 185 186 ierr = PetscPrintf(PETSC_COMM_WORLD,"Initial Temperature = %g\n",(double)sfctemp);CHKERRQ(ierr); /* input surface temperature */ 187 188 deep_grnd_temp = sfctemp - 10; /* set underlying ground layer temperature */ 189 emma = emission(pwat); /* accounts for radiative effects of water vapor */ 190 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); 197 198 /* Converts from Celsius to Kelvin */ 199 sfctemp += 273; 200 airtemp += 273; 201 dewtemp += 273; 202 cloudTemp += 273; 203 deep_grnd_temp += 273; 204 205 /* Calculates initial relative humidity */ 206 x = calcmixingr(dewtemp,pressure1); 207 mixratio = calcmixingr(sfctemp,pressure1); 208 rh = (x/mixratio)*100; 209 210 ierr = PetscPrintf(PETSC_COMM_WORLD,"Initial RH = %.1f percent\n\n",(double)rh);CHKERRQ(ierr); /* prints initial relative humidity */ 211 212 time = 3600*put.time; /* sets amount of timesteps to run model */ 213 214 /* Configure PETSc TS solver */ 215 /*------------------------------------------*/ 216 217 /* Create grid */ 218 ierr = DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,DM_BOUNDARY_PERIODIC,DMDA_STENCIL_STAR,20,20,PETSC_DECIDE,PETSC_DECIDE,dof,1,NULL,NULL,&da);CHKERRQ(ierr); 219 ierr = DMSetFromOptions(da);CHKERRQ(ierr); 220 ierr = DMSetUp(da);CHKERRQ(ierr); 221 ierr = DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);CHKERRQ(ierr); 222 223 /* Define output window for each variable of interest */ 224 ierr = DMDASetFieldName(da,0,"Ts");CHKERRQ(ierr); 225 ierr = DMDASetFieldName(da,1,"Ta");CHKERRQ(ierr); 226 ierr = DMDASetFieldName(da,2,"u");CHKERRQ(ierr); 227 ierr = DMDASetFieldName(da,3,"v");CHKERRQ(ierr); 228 ierr = DMDASetFieldName(da,4,"p");CHKERRQ(ierr); 229 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 */ 245 246 /* set values for MonitorCtx */ 247 usermonitor.drawcontours = PETSC_FALSE; 248 ierr = PetscOptionsHasName(NULL,NULL,"-drawcontours",&usermonitor.drawcontours);CHKERRQ(ierr); 249 if (usermonitor.drawcontours) { 250 PetscReal bounds[] = {1000.0,-1000., -1000.,-1000., 1000.,-1000., 1000.,-1000., 1000,-1000, 100700,100800}; 251 ierr = PetscViewerDrawOpen(PETSC_COMM_WORLD,0,0,0,0,300,300,&usermonitor.drawviewer);CHKERRQ(ierr); 252 ierr = PetscViewerDrawSetBounds(usermonitor.drawviewer,dof,bounds);CHKERRQ(ierr); 253 } 254 usermonitor.interval = 1; 255 ierr = PetscOptionsGetInt(NULL,NULL,"-monitor_interval",&usermonitor.interval,NULL);CHKERRQ(ierr); 256 257 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 258 Extract global vectors from DA; 259 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 260 ierr = DMCreateGlobalVector(da,&T);CHKERRQ(ierr); 261 ierr = VecDuplicate(T,&rhs);CHKERRQ(ierr); /* r: vector to put the computed right hand side */ 262 263 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); 264 ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr); 265 ierr = TSSetType(ts,TSBEULER);CHKERRQ(ierr); 266 ierr = TSSetRHSFunction(ts,rhs,RhsFunc,&user);CHKERRQ(ierr); 267 268 /* Set Jacobian evaluation routine - use coloring to compute finite difference Jacobian efficiently */ 269 ierr = DMSetMatType(da,MATAIJ);CHKERRQ(ierr); 270 ierr = DMCreateMatrix(da,&J);CHKERRQ(ierr); 271 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 272 if (use_coloring) { 273 ISColoring iscoloring; 274 ierr = DMCreateColoring(da,IS_COLORING_GLOBAL,&iscoloring);CHKERRQ(ierr); 275 ierr = MatFDColoringCreate(J,iscoloring,&matfdcoloring);CHKERRQ(ierr); 276 ierr = MatFDColoringSetFromOptions(matfdcoloring);CHKERRQ(ierr); 277 ierr = MatFDColoringSetUp(J,iscoloring,matfdcoloring);CHKERRQ(ierr); 278 ierr = ISColoringDestroy(&iscoloring);CHKERRQ(ierr); 279 ierr = MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))SNESTSFormFunction,ts);CHKERRQ(ierr); 280 ierr = SNESSetJacobian(snes,J,J,SNESComputeJacobianDefaultColor,matfdcoloring);CHKERRQ(ierr); 281 } else { 282 ierr = SNESSetJacobian(snes,J,J,SNESComputeJacobianDefault,NULL);CHKERRQ(ierr); 283 } 284 285 /* Define what to print for ts_monitor option */ 286 ierr = PetscOptionsHasName(NULL,NULL,"-monitor_off",&monitor_off);CHKERRQ(ierr); 287 if (!monitor_off) { 288 ierr = TSMonitorSet(ts,Monitor,&usermonitor,NULL);CHKERRQ(ierr); 289 } 290 ierr = FormInitialSolution(da,T,&user);CHKERRQ(ierr); 291 dt = TIMESTEP; /* initial time step */ 292 ftime = TIMESTEP*time; 293 ierr = PetscPrintf(PETSC_COMM_WORLD,"time %D, ftime %g hour, TIMESTEP %g\n",time,(double)(ftime/3600),(double)dt);CHKERRQ(ierr); 294 295 ierr = TSSetTimeStep(ts,dt);CHKERRQ(ierr); 296 ierr = TSSetMaxSteps(ts,time);CHKERRQ(ierr); 297 ierr = TSSetMaxTime(ts,ftime);CHKERRQ(ierr); 298 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr); 299 ierr = TSSetSolution(ts,T);CHKERRQ(ierr); 300 ierr = TSSetDM(ts,da);CHKERRQ(ierr); 301 302 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 303 Set runtime options 304 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 305 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 306 307 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 308 Solve nonlinear system 309 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 310 ierr = TSSolve(ts,T);CHKERRQ(ierr); 311 ierr = TSGetSolveTime(ts,&ftime);CHKERRQ(ierr); 312 ierr = TSGetStepNumber(ts,&steps);CHKERRQ(ierr); 313 ierr = PetscPrintf(PETSC_COMM_WORLD,"Solution T after %g hours %D steps\n",(double)(ftime/3600),steps);CHKERRQ(ierr); 314 315 if (matfdcoloring) {ierr = MatFDColoringDestroy(&matfdcoloring);CHKERRQ(ierr);} 316 if (usermonitor.drawcontours) { 317 ierr = PetscViewerDestroy(&usermonitor.drawviewer);CHKERRQ(ierr); 318 } 319 ierr = MatDestroy(&J);CHKERRQ(ierr); 320 ierr = VecDestroy(&T);CHKERRQ(ierr); 321 ierr = VecDestroy(&rhs);CHKERRQ(ierr); 322 ierr = TSDestroy(&ts);CHKERRQ(ierr); 323 ierr = DMDestroy(&da);CHKERRQ(ierr); 324 325 ierr = 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 { 334 PetscFunctionBeginUser; 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 PetscFunctionReturn(0); 337 } 338 339 PetscErrorCode calcfluxa(PetscScalar sfctemp, PetscScalar airtemp, PetscScalar emma, PetscScalar *flux) /* this function is not currently called upon */ 340 { 341 PetscScalar emm = 0.001; 342 343 PetscFunctionBeginUser; 344 *flux = SIG*(-emm*(PetscPowScalarInt(airtemp,4))); /* calculates flux usinge Stefan-Boltzmann relation */ 345 PetscFunctionReturn(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 */ 352 353 PetscFunctionBeginUser; 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 PetscFunctionReturn(0); 357 } 358 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 */ 370 371 PetscFunctionBeginUser; 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 */ 379 380 *latentheat = density*wndmix*beta*lhcnst*(q - qs); /* calculates latent heat flux */ 381 PetscFunctionReturn(0); 382 } 383 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 */ 390 391 PetscFunctionBeginUser; 392 /* mixratio = calcmixingr(sfctemp,pressure1); */ 393 394 /* initialize poisson constant */ 395 kdry = 0.2854; 396 /* kmoist = 0.2854*(1 - 0.24*mixratio); */ 397 398 pavg = ((0.7*pressure1)+pressure2)/2; /* calculates simple average press */ 399 *pottemp = temp*(PetscPowScalar((pressure1/pavg),kdry)); /* calculates potential temperature */ 400 PetscFunctionReturn(0); 401 } 402 extern PetscScalar calcmixingr(PetscScalar dtemp, PetscScalar pressure1) 403 { 404 PetscScalar e; /* vapor pressure */ 405 PetscScalar mixratio; /* mixing ratio */ 406 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 */ 412 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 } 421 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 */ 428 429 PetscFunctionBeginUser; 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 PetscFunctionReturn(0); 433 } 434 extern PetscScalar emission(PetscScalar pwat) 435 { 436 PetscScalar emma; 437 438 emma = 0.725 + 0.17*PetscLog10Real(PetscRealPart(pwat)); 439 440 return emma; 441 } 442 extern PetscScalar cloud(PetscScalar fract) 443 { 444 PetscScalar emma = 0; 445 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; 483 484 PetscFunctionBegin; 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; 490 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; 494 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; 498 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; 502 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; 506 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; 510 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; 514 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; 518 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; 522 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 PetscFunctionReturn(0); 527 } 528 529 /* ------------------------------------------------------------------- */ 530 PetscErrorCode FormInitialSolution(DM da,Vec Xglobal,void *ctx) 531 { 532 PetscErrorCode ierr; 533 AppCtx *user = (AppCtx*)ctx; /* user-defined application context */ 534 PetscInt i,j,xs,ys,xm,ym,Mx,My; 535 Field **X; 536 537 PetscFunctionBeginUser; 538 ierr = 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);CHKERRQ(ierr); 540 541 /* Get pointers to vector data */ 542 ierr = DMDAVecGetArray(da,Xglobal,&X);CHKERRQ(ierr); 543 544 /* Get local grid boundaries */ 545 ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr); 546 547 /* Compute function over the locally owned part of the grid */ 548 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 } 572 573 /* Restore vectors */ 574 ierr = DMDAVecRestoreArray(da,Xglobal,&X);CHKERRQ(ierr); 575 PetscFunctionReturn(0); 576 } 577 578 /* 579 RhsFunc - Evaluates nonlinear function F(u). 580 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() 587 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; 595 PetscErrorCode ierr; 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; 617 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; 623 624 PetscFunctionBeginUser; 625 ierr = DMGetLocalVector(da,&localT);CHKERRQ(ierr); 626 ierr = 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);CHKERRQ(ierr); 627 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; */ 630 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 ierr = DMGlobalToLocalBegin(da,Xglobal,INSERT_VALUES,localT);CHKERRQ(ierr); 638 ierr = DMGlobalToLocalEnd(da,Xglobal,INSERT_VALUES,localT);CHKERRQ(ierr); 639 640 /* Get pointers to vector data */ 641 ierr = DMDAVecGetArrayRead(da,localT,&X);CHKERRQ(ierr); 642 ierr = DMDAVecGetArray(da,F,&Frhs);CHKERRQ(ierr); 643 644 /* Get local grid boundaries */ 645 ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr); 646 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; */ 653 654 sfctemp1 = (double)Ts; 655 ierr = calcfluxs(sfctemp1,airtemp,emma,fract,Tc,&fsfc1);CHKERRQ(ierr); /* calculates surface net radiative flux */ 656 ierr = sensibleflux(sfctemp1,airtemp,wind,&sheat);CHKERRQ(ierr); /* calculate sensible heat flux */ 657 ierr = latentflux(sfctemp1,dewtemp,wind,pressure1,&latentheat);CHKERRQ(ierr); /* calculates latent heat flux */ 658 ierr = calc_gflux(sfctemp1,deep_grnd_temp,&groundflux);CHKERRQ(ierr); /* calculates flux from earth below surface soil layer by conduction */ 659 ierr = calcfluxa(sfctemp1,airtemp,emma,&Ra);CHKERRQ(ierr); /* 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 */ 661 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 */ 666 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 */ 670 671 /* Solve governing equations */ 672 /* P = p*Rd*Ts; */ 673 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; 680 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; 687 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); 694 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) */ 698 699 Frhs[j][i].Ta = Ra/Cp; /* dTa/dt time change of air temperature */ 700 } 701 } 702 703 /* Restore vectors */ 704 ierr = DMDAVecRestoreArrayRead(da,localT,&X);CHKERRQ(ierr); 705 ierr = DMDAVecRestoreArray(da,F,&Frhs);CHKERRQ(ierr); 706 ierr = DMRestoreLocalVector(da,&localT);CHKERRQ(ierr); 707 PetscFunctionReturn(0); 708 } 709 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; 717 718 PetscFunctionBeginUser; 719 ierr = VecNorm(T,NORM_INFINITY,&norm);CHKERRQ(ierr); 720 721 if (step%user->interval == 0) { 722 ierr = VecGetArrayRead(T,&array);CHKERRQ(ierr); 723 ierr = 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 ierr = VecRestoreArrayRead(T,&array);CHKERRQ(ierr); 725 } 726 727 if (user->drawcontours) { 728 ierr = VecView(T,viewer);CHKERRQ(ierr); 729 } 730 PetscFunctionReturn(0); 731 } 732 733 /*TEST 734 735 build: 736 requires: !complex !single 737 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 743 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 751 752 TEST*/ 753