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