258
QH2=0;// reduced thermal flux
259
GRADTHETA1=0;
260
GRADTHETA2=0;// reduced thermal gradient
261
GRADPHI1=0;
262
GRADPHI2=0;// reduced gradient of PHI from Clapeyron
263
Z=0;// reduced thickness
264
ZN=0;//reduced location of new lens
265
ZI=0;//reduced location of old lens
266
DZ=0;// reduced layer thickness
267
UW=0;//reduced water pressure
268
DUW=0;// reduced change in water pressure
269
// across layer
270
UI=PHI;// reduced ice pressure
271
CHI=0;//reduced value of the stress partition factor
272
UN=0;// reduced neutral stress
273
UNold=0;// previous value of reduced neutral stress
274
UNmax=-1000;// maximum value of neutral stress
275
THETA=0;
276
DTHETA=0;
277
GRADPROB=0;
278
LAYERPROB=0;
279
InitStress();
280
InitTol();
281
}
282
283
284
void Trigidice::LAYER(double rphi, double rdphi){
285
DW = W(DSAT(rphi))-W(DSAT(rphi+rdphi));
286
QW2 = QW1 + (VB-Y*(VB+VI))*DW;
287
GRADUW2 = FW - QW2/RKW(rphi+rdphi);
288
QH2 = QH1 + H*(DW*VB + (QW2 - QW1));
289
GRADTHETA2 = -QH2/RKH(W(DSAT(rphi+rdphi)));
290
GRADPHI2 = (Y-1)* GRADUW2 - Y*H* GRADTHETA2;
291
if ((GRADPHI1<0)||(GRADPHI2< 0)) {
292
GRADPROB=1;
293
DZ=1;
294
} else
295
DZ = rdphi/(sqrt(GRADPHI1*GRADPHI2));
296
DUW = FW*DZ - DZ*(QW1/RKW(rphi) +
297
QW2/(RKW(rphi+rdphi)))/2;
298
DTHETA = ((Y-1)*DUW-rdphi)/(Y*H);
299
// reset variables for beginning of next layer
300
QW1=QW2;
301
GRADUW1=GRADUW2;
302
QH1=QH2;
303
GRADTHETA1=GRADTHETA2;
304
GRADPHI1=GRADPHI2;
305
}
306
307
void Trigidice::PROFILE(){
308
LAYER(PHI,DPHI);
309
PHI=PHI+DPHI;
310
DPHI = resolution/DZ;
311
THETA=THETA+DTHETA;
23