OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
mat104c_nodam_nice.F File Reference
#include "implicit_f.inc"
#include "vectorize.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine mat104c_nodam_nice (nel, ngl, nuparam, nuvar, time, timestep, uparam, uvar, jthe, off, gs, rho, pla, dpla, epsd, soundsp, depsxx, depsyy, depsxy, depsyz, depszx, sigoxx, sigoyy, sigoxy, sigoyz, sigozx, signxx, signyy, signxy, signyz, signzx, thkly, thk, sigy, et, tempel, dplanl, temp, seq, inloc)

Function/Subroutine Documentation

◆ mat104c_nodam_nice()

subroutine mat104c_nodam_nice ( integer nel,
integer, dimension(nel), intent(in) ngl,
integer nuparam,
integer nuvar,
time,
timestep,
intent(in) uparam,
intent(inout) uvar,
integer jthe,
intent(inout) off,
intent(in) gs,
intent(in) rho,
intent(inout) pla,
intent(inout) dpla,
intent(inout) epsd,
intent(out) soundsp,
intent(in) depsxx,
intent(in) depsyy,
intent(in) depsxy,
intent(in) depsyz,
intent(in) depszx,
intent(in) sigoxx,
intent(in) sigoyy,
intent(in) sigoxy,
intent(in) sigoyz,
intent(in) sigozx,
intent(out) signxx,
intent(out) signyy,
intent(out) signxy,
intent(out) signyz,
intent(out) signzx,
intent(in) thkly,
intent(inout) thk,
intent(out) sigy,
intent(out) et,
intent(in) tempel,
intent(in) dplanl,
intent(inout) temp,
intent(inout) seq,
integer inloc )

Definition at line 28 of file mat104c_nodam_nice.F.

37 !=======================================================================
38 ! Implicit types
39 !=======================================================================
40#include "implicit_f.inc"
41 !=======================================================================
42 ! Common
43 !=======================================================================
44 !=======================================================================
45 ! Dummy arguments
46 !=======================================================================
47 INTEGER NEL,NUPARAM,NUVAR,JTHE,INLOC
48 INTEGER ,DIMENSION(NEL), INTENT(IN) :: NGL
49 my_real
50 . time,timestep
51 my_real,DIMENSION(NUPARAM), INTENT(IN) ::
52 . uparam
53 my_real,DIMENSION(NEL), INTENT(IN) ::
54 . rho,tempel,dplanl,
55 . depsxx,depsyy,depsxy,depsyz,depszx,
56 . sigoxx,sigoyy,sigoxy,sigoyz,sigozx,
57 . gs , thkly
58c
59 my_real ,DIMENSION(NEL), INTENT(OUT) ::
60 . soundsp,sigy,et,
61 . signxx,signyy,signxy,signyz,signzx
62c
63 my_real ,DIMENSION(NEL), INTENT(INOUT) ::
64 . pla,dpla,epsd,off,thk,temp,seq
65 my_real ,DIMENSION(NEL,NUVAR), INTENT(INOUT) ::
66 . uvar
67c
68 !=======================================================================
69 ! Local Variables
70 !=======================================================================
71 INTEGER I,K,II,NINDX,INDEX(NEL),NICE
72c
73 my_real
74 . young,bulk,lam,g,g2,nu,cdr,kdr,hard,yld0,qvoce,bvoce,jcc,
75 . epsp0,mtemp,tini,tref,eta,cp,dpis,dpad,asrate,afiltr,a11,a12,
76 . nnu, normsig
77 my_real
78 . h,ldav,trdfds,sigvm,sigdr2,yld2i,omega,
79 . dtemp,fcosh,fsinh,dpdt,dlam,ddep,depxx,depyy
80 my_real
81 . dpxx,dpyy,dpxy,dpzz,dsdrdj2,dsdrdj3,
82 . dj3dsxx,dj3dsyy,dj3dsxy,dj3dszz,
83 . dj2dsxx,dj2dsyy,dj2dsxy,
84 . dfdsxx,dfdsyy,dfdsxy,
85 . normxx,normyy,normxy,normzz,
86 . denom,dphi,dphi_dtrsig,sig_dfdsig,dfdsig2,
87 . dphi_dsig,dphi_dyld,dphi_dfdr,dpdt_nl,
88 . dyld_dpla,dyld_dtemp,dtemp_dlam,dlam_nl
89c
90 my_real, DIMENSION(NEL) ::
91 . dsigxx,dsigyy,dsigxy,trsig,trdep,
92 . sxx,syy,sxy,szz,sigm,j2,j3,sigdr,yld,weitemp,
93 . hardp,fhard,frate,ftherm,fdr,phi0,phi,dpla_dlam,
94 . dezz
95 !=======================================================================
96 ! DRUCKER - VOCE - JOHNSON-COOK MATERIAL LAW
97 !=======================================================================
98 !UVAR(1) PHI YIELD FUNCTION VALUE
99 !=======================================================================
100c
101 !=======================================================================
102 ! - INITIALISATION OF COMPUTATION ON TIME STEP
103 !=======================================================================
104 ! Recovering model parameters
105 ! Elastic parameters
106 young = uparam(1) ! Young modulus
107 bulk = uparam(2) ! Bulk modulus
108 g = uparam(3) ! Shear modulus
109 g2 = uparam(4) ! 2*Shear modulus
110 lam = uparam(5) ! Lambda Hooke parameter
111 nu = uparam(6) ! Poisson ration
112 nnu = uparam(7) ! NU/(1-NU)
113 a11 = uparam(9) ! Diagonal term, elastic matrix in plane stress
114 a12 = uparam(10) ! Non-diagonal term, elastic matrix in plane stress
115 ! plastic criterion and hardening parameters [drucker, 1948]
116 nice = nint(uparam(11)) ! Choice of the Nice method
117 cdr = uparam(12) ! Drucker coefficient
118 kdr = uparam(13) ! Drucker 1/K coefficient
119 tini = uparam(14) ! Initial temperature
120 hard = uparam(15) ! Linear hardening
121 yld0 = uparam(16) ! Initial yield stress
122 qvoce = uparam(17) ! 1st Voce parameter
123 bvoce = uparam(18) ! 2nd Voce parameter
124 ! Strain-rate dependence parameters
125 jcc = uparam(20) ! Johnson-Cook strain rate coefficient
126 epsp0 = uparam(21) ! Johnson-Cook reference strain rate
127 ! Thermal softening and self-heating parameters
128 mtemp = uparam(22) ! Thermal softening slope
129 tref = uparam(23) ! Reference temperature
130 eta = uparam(24) ! Taylor-Quinney coefficient
131 cp = uparam(25) ! Thermal mass capacity
132 dpis = uparam(26) ! Isothermal plastic strain rate
133 dpad = uparam(27) ! Adiabatic plastic strain rate
134 ! Plastic strain-rate filtering parameters
135 asrate = uparam(28) ! Plastic strain rate filtering frequency
136 afiltr = asrate*timestep/(asrate*timestep + one)
137c
138 ! Recovering internal variables
139 DO i=1,nel
140 IF (off(i) < em01) off(i) = zero
141 IF (off(i) < one) off(i) = off(i)*four_over_5
142 ! User inputs
143 phi0(i) = uvar(i,1) ! Previous yield function value
144 ! Standard inputs
145 dpla(i) = zero ! Initialization of the plastic strain increment
146 et(i) = one ! Initialization of hourglass coefficient
147 hardp(i) = zero ! Initialization of hardening modulus
148 ENDDO
149c
150 ! Initialization of temperature and self-heating weight factor
151 IF (time == zero) THEN
152 temp(1:nel) = tini
153 ENDIF
154 IF (cp > zero) THEN
155 IF (jthe == 0) THEN
156 IF (inloc == 0) THEN
157 DO i=1,nel
158 ! Computation of the weight factor
159 IF (epsd(i) < dpis) THEN
160 weitemp(i) = zero
161 ELSEIF (epsd(i) > dpad) THEN
162 weitemp(i) = one
163 ELSE
164 weitemp(i) = ((epsd(i)-dpis)**2 )
165 . * (three*dpad - two*epsd(i) - dpis)
166 . / ((dpad-dpis)**3)
167 ENDIF
168 ENDDO
169 ENDIF
170 ELSE
171 temp(1:nel) = tempel(1:nel)
172 ENDIF
173 ENDIF
174c
175 ! Computation of the initial yield stress
176 DO i = 1,nel
177 ! a) - Hardening law
178 fhard(i) = yld0 + hard*pla(i) + qvoce*(one-exp(-bvoce*pla(i)))
179 ! b) - Correction factor for strain-rate dependence (Johnson-Cook)
180 frate(i) = one
181 IF (epsd(i) > epsp0) frate(i) = frate(i) + jcc*log(epsd(i)/epsp0)
182 ! c) - Correction factor for thermal softening
183 ftherm(i) = one
184 IF (cp > zero) ftherm(i) = one - mtemp * (temp(i) - tref)
185 ! d) - Computation of the yield function and value check
186 yld(i) = fhard(i)*frate(i)*ftherm(i)
187 ! e) - Checking values
188 yld(i) = max(em10, yld(i))
189 ENDDO
190c
191 !========================================================================
192 ! - COMPUTATION OF TRIAL VALUES
193 !========================================================================
194 DO i=1,nel
195c
196 ! Computation of the trial stress tensor
197 signxx(i) = sigoxx(i) + a11*depsxx(i) + a12*depsyy(i)
198 signyy(i) = sigoyy(i) + a11*depsyy(i) + a12*depsxx(i)
199 signxy(i) = sigoxy(i) + depsxy(i)*g
200 signyz(i) = sigoyz(i) + depsyz(i)*gs(i)
201 signzx(i) = sigozx(i) + depszx(i)*gs(i)
202 ! Computation of the trace of the trial stress tensor
203 trsig(i) = signxx(i) + signyy(i)
204 sigm(i) = -trsig(i) * third
205 ! Computation of the deviatoric trial stress tensor
206 sxx(i) = signxx(i) + sigm(i)
207 syy(i) = signyy(i) + sigm(i)
208 szz(i) = sigm(i)
209 sxy(i) = signxy(i)
210 dezz(i) = -nnu*depsxx(i) - nnu*depsyy(i)
211 ! Second deviatoric invariant
212 j2(i) = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2 ) + sxy(i)**2
213 ! Third deviatoric invariant
214 j3(i) = sxx(i)*syy(i)*szz(i) - szz(i)*sxy(i)**2
215 ! Drucker equivalent stress
216 fdr(i) = j2(i)**3 - cdr*(j3(i)**2)
217 ! Checking equivalent stress values
218 IF (fdr(i) > zero) THEN
219 sigdr(i) = kdr * exp((one/six)*log(fdr(i))) ! FDR(I)**(1/6)
220 ELSE
221 sigdr(i) = zero
222 ENDIF
223 ENDDO
224c
225 !========================================================================
226 ! - COMPUTATION OF YIELD FONCTION
227 !========================================================================
228 phi(1:nel) = (sigdr(1:nel) / yld(1:nel))**2 - one
229c
230 ! Checking plastic behavior for all elements
231 nindx = 0
232 DO i=1,nel
233 IF (phi(i) >= zero .AND. off(i) == one) THEN
234 nindx=nindx+1
235 index(nindx)=i
236 ENDIF
237 ENDDO
238c
239 !====================================================================
240 ! - PLASTIC CORRECTION WITH NICE - EXPLICIT METHOD
241 !====================================================================
242#include "vectorize.inc"
243 DO ii=1,nindx
244c
245 ! Number of the element with plastic behaviour
246 i = index(ii)
247c
248 ! Computation of the trial stress increment
249 dsigxx(i) = a11*depsxx(i) + a12*depsyy(i)
250 dsigyy(i) = a11*depsyy(i) + a12*depsxx(i)
251 dsigxy(i) = depsxy(i)*g
252 dlam = zero
253c
254 ! Computation of Drucker equivalent stress of the previous stress tensor
255 trsig(i) = sigoxx(i) + sigoyy(i)
256 sigm(i) = -trsig(i) * third
257 sxx(i) = sigoxx(i) + sigm(i)
258 syy(i) = sigoyy(i) + sigm(i)
259 szz(i) = sigm(i)
260 sxy(i) = sigoxy(i)
261 j2(i) = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2) + sxy(i)**2
262 j3(i) = sxx(i)*syy(i)*szz(i) - szz(i)*sxy(i)**2
263 fdr(i) = j2(i)**3 - cdr*(j3(i)**2)
264 sigdr(i) = kdr * exp((one/six)*log(fdr(i)))
265c
266 ! Note : in this part, the purpose is to compute in one iteration
267 ! a plastic multiplier allowing to update internal variables to satisfy
268 ! the consistency condition.
269 ! Its expression is : DLAMBDA = - (PHI_OLD + DPHI) / DPHI_DLAMBDA
270 ! -> PHI_OLD : old value of yield function (known)
271 ! -> DPHI : yield function prediction (to compute)
272 ! -> DPHI_DLAMBDA : derivative of PHI with respect to DLAMBDA by taking
273 ! into account of internal variables kinetic :
274 ! plasticity, temperature and damage (to compute)
275c
276 ! 1 - Computation of DPHI_DSIG the normal to the yield surface
277 !-------------------------------------------------------------
278c
279 ! Derivative with respect to the equivalent stress
280 yld2i = one/(yld(i)**2)
281 dphi_dsig = two*sigdr(i)*yld2i
282c
283 ! Computation of the Eulerian norm of the stress tensor
284 normsig = sqrt(sigoxx(i)*sigoxx(i)
285 . + sigoyy(i)*sigoyy(i)
286 . + two*sigoxy(i)*sigoxy(i))
287 normsig = max(normsig,one)
288c
289 ! Derivative with respect to Fdr
290 fdr(i) = (j2(i)/(normsig**2))**3 - cdr*((j3(i)/(normsig**3))**2)
291 dphi_dfdr = dphi_dsig*kdr*(one/six)*exp(-(five/six)*log(fdr(i)))
292 dsdrdj2 = dphi_dfdr*three*(j2(i)/(normsig**2))**2
293 dsdrdj3 = -dphi_dfdr*two*cdr*(j3(i)/(normsig**3))
294 ! dJ3/dSig
295 dj3dsxx = two_third*(syy(i)*szz(i))/(normsig**2)
296 . - third*(sxx(i)*szz(i))/(normsig**2)
297 . - third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
298 dj3dsyy = - third*(syy(i)*szz(i))/(normsig**2)
299 . + two_third*(sxx(i)*szz(i))/(normsig**2)
300 . - third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
301 dj3dszz = - third*(syy(i)*szz(i))/(normsig**2)
302 . - third*(sxx(i)*szz(i))/(normsig**2)
303 . + two_third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
304 dj3dsxy = two*(sxx(i)*sxy(i) + sxy(i)*syy(i))/(normsig**2)
305 ! dPhi/dSig
306 normxx = dsdrdj2*sxx(i)/normsig + dsdrdj3*dj3dsxx
307 normyy = dsdrdj2*syy(i)/normsig + dsdrdj3*dj3dsyy
308 normzz = dsdrdj2*szz(i)/normsig + dsdrdj3*dj3dszz
309 normxy = two*dsdrdj2*sxy(i)/normsig + dsdrdj3*dj3dsxy
310c
311 ! Restoring previous value of the yield function
312 phi(i) = phi0(i)
313c
314 ! Computation of yield surface trial increment DPHI
315 dphi = normxx * dsigxx(i)
316 . + normyy * dsigyy(i)
317 . + normxy * dsigxy(i)
318c
319 ! 2 - computation of dphi_dlambda
320 !---------------------------------------------------------
321c
322 ! a) Derivative with respect stress increments tensor DSIG
323 ! --------------------------------------------------------
324 dfdsig2 = normxx * (a11*normxx + a12*normyy)
325 . + normyy * (a11*normyy + a12*normxx)
326 . + normxy * normxy * g
327c
328 ! b) Derivatives with respect to plastic strain P
329 ! ------------------------------------------------
330c
331 ! i) Derivative of the yield stress with respect to plastic strain dYLD / dPLA
332 ! ----------------------------------------------------------------------------
333 hardp(i) = hard + qvoce*bvoce*exp(-bvoce*pla(i))
334 dyld_dpla = hardp(i)*frate(i)*ftherm(i)
335c
336 ! ii) Derivative of dPLA with respect to DLAM
337 ! -------------------------------------------
338 sig_dfdsig = sigoxx(i) * normxx
339 . + sigoyy(i) * normyy
340 . + sigoxy(i) * normxy
341 dpla_dlam(i) = sig_dfdsig / yld(i)
342c
343 ! c) Derivatives with respect to the temperature TEMP
344 ! ---------------------------------------------------
345 IF (jthe == 0 .AND. cp > zero .AND. inloc == 0) THEN
346 ! i) Derivative of the yield stress with respect to temperature dYLD / dTEMP
347 ! ---------------------------------------------------------------------------
348 dyld_dtemp = -fhard(i)*frate(i)*mtemp
349 ! ii) Derivative of the temperature TEMP with respect to DLAM
350 ! -----------------------------------------------------------
351 dtemp_dlam = weitemp(i)*(eta/(rho(i)*cp))*sig_dfdsig
352 ELSE
353 dyld_dtemp = zero
354 dtemp_dlam = zero
355 ENDIF
356c
357 ! d) Derivative with respect to the yield stress
358 ! ----------------------------------------------
359 sigdr2 = sigdr(i)**2
360 dphi_dyld = -two*sigdr2/(yld(i)**3)
361c
362 ! e) Derivative of PHI with respect to DLAM ( = -DENOM)
363 denom = dfdsig2 - (dphi_dyld * dyld_dpla * dpla_dlam(i))
364 IF (jthe == 0 .AND. cp > zero .AND. inloc == 0) THEN
365 denom = denom - (dphi_dyld * dyld_dtemp * dtemp_dlam)
366 ENDIF
367 denom = sign(max(abs(denom),em20) ,denom)
368c
369 ! 3 - Computation of plastic multiplier and variables update
370 !----------------------------------------------------------
371c
372 ! Computation of the plastic multiplier
373 dlam = (dphi + phi(i)) / denom
374 dlam = max(dlam, zero)
375c
376 ! Plastic strains tensor update
377 dpxx = dlam * normxx
378 dpyy = dlam * normyy
379 dpzz = dlam * normzz
380 dpxy = dlam * normxy
381c
382 ! Elasto-plastic stresses update
383 signxx(i) = sigoxx(i) + dsigxx(i) - (a11*dpxx + a12*dpyy)
384 signyy(i) = sigoyy(i) + dsigyy(i) - (a11*dpyy + a12*dpxx)
385 signxy(i) = sigoxy(i) + dsigxy(i) - dpxy*g
386 trsig(i) = signxx(i) + signyy(i)
387 sigm(i) = -trsig(i) * third
388 sxx(i) = signxx(i) + sigm(i)
389 syy(i) = signyy(i) + sigm(i)
390 szz(i) = sigm(i)
391 sxy(i) = signxy(i)
392c
393 ! Cumulated plastic strain and strain rate update
394 ddep = (dlam/yld(i))*sig_dfdsig
395 dpla(i) = dpla(i) + ddep
396 pla(i) = pla(i) + ddep
397c
398 ! Drucker equivalent stress update
399 j2(i) = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2 ) + sxy(i)**2
400 j3(i) = sxx(i) * syy(i) * szz(i) - szz(i) * sxy(i)**2
401 fdr(i) = j2(i)**3 - cdr*(j3(i)**2)
402 sigdr(i) = kdr * exp((one/six)*log(fdr(i)))
403c
404 ! Temperature update
405 IF (jthe == 0 .AND. cp > zero .AND. inloc == 0) THEN
406 dtemp = weitemp(i)*yld(i)*ddep*eta/(rho(i)*cp)
407 temp(i) = temp(i) + dtemp
408 ftherm(i) = one - mtemp* (temp(i) - tref)
409 ENDIF
410c
411 ! Hardening law update
412 fhard(i) = yld0 + hard * pla(i) + qvoce*(one-exp(-bvoce*pla(i)))
413c
414 ! Yield stress update
415 yld(i) = fhard(i) * frate(i) * ftherm(i)
416 yld(i) = max(yld(i), em10)
417c
418 ! Yield function value update
419 sigdr2 = sigdr(i)**2
420 yld2i = one / yld(i)**2
421 phi(i) = sigdr2 * yld2i - one
422c
423 ! Transverse strain update
424 IF (inloc == 0) THEN
425 dezz(i) = dezz(i) + nnu*dpxx + nnu*dpyy + dpzz
426 ENDIF
427 ENDDO
428 !===================================================================
429 ! - END OF PLASTIC CORRECTION WITH NICE - EXPLICIT METHOD
430 !===================================================================
431c
432 ! Storing new values
433 DO i=1,nel
434 ! USR Outputs & coefficient for hourglass
435 IF (dpla(i) > zero) THEN
436 uvar(i,1) = phi(i) ! Yield function value
437 et(i) = hardp(i)*frate(i) / (hardp(i)*frate(i) + young)
438 ELSE
439 uvar(i,1) = zero
440 et(i) = one
441 ENDIF
442 seq(i) = sigdr(i)
443 ! Plastic strain-rate (filtered)
444 dpdt = dpla(i) / max(em20,timestep)
445 epsd(i) = afiltr * dpdt + (one - afiltr) * epsd(i)
446 ! Computation of the sound speed
447 soundsp(i) = sqrt(a11/rho(i))
448 ! Storing the yield stress
449 sigy(i) = yld(i)
450 ! Non-local thickness variation + temperature
451 IF ((inloc > 0).AND.(off(i) == one).AND.(dplanl(i)>zero)) THEN
452 yld2i = one/(yld(i)**2)
453 dphi_dsig = two*sigdr(i)*yld2i
454 normsig = sqrt(signxx(i)*signxx(i)
455 . + signyy(i)*signyy(i)
456 . + two*signxy(i)*signxy(i))
457 normsig = max(normsig,one)
458 fdr(i) = (j2(i)/(normsig**2))**3 - cdr*((j3(i)/(normsig**3))**2)
459 IF (fdr(i) > zero) THEN
460 dphi_dfdr = dphi_dsig*kdr*(one/six)*exp(-(five/six)*log(fdr(i)))
461 ELSE
462 dphi_dfdr = zero
463 ENDIF
464 dsdrdj2 = dphi_dfdr*three*(j2(i)/(normsig**2))**2
465 dsdrdj3 = -dphi_dfdr*two*cdr*(j3(i)/(normsig**3))
466 dj3dsxx = two_third*(syy(i)*szz(i))/(normsig**2)
467 . - third*(sxx(i)*szz(i))/(normsig**2)
468 . - third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
469 dj3dsyy = - third*(syy(i)*szz(i))/(normsig**2)
470 . + two_third*(sxx(i)*szz(i))/(normsig**2)
471 . - third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
472 dj3dszz = - third*(syy(i)*szz(i))/(normsig**2)
473 . - third*(sxx(i)*szz(i))/(normsig**2)
474 . + two_third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
475 dj3dsxy = two*(sxx(i)*sxy(i) + sxy(i)*syy(i))/(normsig**2)
476 normxx = dsdrdj2*sxx(i)/normsig + dsdrdj3*dj3dsxx
477 normyy = dsdrdj2*syy(i)/normsig + dsdrdj3*dj3dsyy
478 normzz = dsdrdj2*szz(i)/normsig + dsdrdj3*dj3dszz
479 normxy = two*dsdrdj2*sxy(i)/normsig + dsdrdj3*dj3dsxy
480 sig_dfdsig = signxx(i) * normxx
481 . + signyy(i) * normyy
482 . + signxy(i) * normxy
483 IF (sig_dfdsig > em01) THEN
484 dlam_nl = (yld(i)*dplanl(i))/sig_dfdsig
485 dezz(i) = dezz(i) + nnu*(dlam_nl*normxx)
486 . + nnu*(dlam_nl*normyy)
487 . + dlam_nl*normzz
488 ENDIF
489 IF (cp > zero .AND. jthe == 0) THEN
490 ! Computation of the weight factor
491 dpdt_nl = dplanl(i)/max(timestep,em20)
492 IF (dpdt_nl < dpis) THEN
493 weitemp(i) = zero
494 ELSEIF (dpdt_nl > dpad) THEN
495 weitemp(i) = one
496 ELSE
497 weitemp(i) = ((dpdt_nl-dpis)**2 )
498 . * (three*dpad - two*dpdt_nl - dpis)
499 . / ((dpad-dpis)**3)
500 ENDIF
501 dtemp = weitemp(i)*dplanl(i)*yld(i)*eta/(rho(i)*cp)
502 temp(i) = temp(i) + dtemp
503 ENDIF
504 ENDIF
505 ! Computation of the thickness variation
506 thk(i) = thk(i) + dezz(i)*thkly(i)*off(i)
507 ENDDO
508c
my_real
#define my_real
Definition cppsort.cpp:32
max
#define max(a, b)
Definition macros.h:21