Functions/Subroutines | |
| subroutine | heatcond1d_vstep (t, qs, tsfc, dx, dt, k, rhocp, bc, bctype, debug) |
| subroutine | heatcond1d (t, qs, dx, dt, k, rhocp, bc, bctype) |
| recursive subroutine | heatcond1d_ext (t, qs, tsfc, dx, dt, k, rhocp, bc, bctype, debug) |
Function/Subroutine Documentation
◆ heatcond1d()
| subroutine heatflux::heatcond1d | ( | real(kind(1d0)), dimension(:), intent(inout) | t, |
| real(kind(1d0)), intent(out) | qs, | ||
| real(kind(1d0)), dimension(:), intent(in) | dx, | ||
| real(kind(1d0)), intent(in) | dt, | ||
| real(kind(1d0)), dimension(:), intent(in) | k, | ||
| real(kind(1d0)), dimension(:), intent(in) | rhocp, | ||
| real(kind(1d0)), dimension(2), intent(in) | bc, | ||
| logical, dimension(2), intent(in) | bctype ) |
Definition at line 480 of file suews_phys_estm.f95.
481 REAL(KIND(1D0)), INTENT(inout) :: T(:)
482 REAL(KIND(1D0)), INTENT(in) :: dx(:), dt, k(:), rhocp(:), bc(2)
483 REAL(KIND(1D0)), INTENT(out) :: Qs
484 LOGICAL, INTENT(in) :: bctype(2)
485 INTEGER :: i, n !,j !!!!!FO!!!!!
486 REAL(KIND(1D0)), ALLOCATABLE :: w(:), a(:), T1(:)
487 n = SIZE(t)
488 ALLOCATE (w(0:n), a(n), t1(n))
489 ! w = interface temperature
490 w(1:n) = t
491 w(0) = bc(1); w(n) = bc(2)
492 !convert from flux to equivalent temperature, not exact
493 ! F = k dT/dX => dx*F/k + T(i+1) = T(i)
494 IF (bctype(1)) w(0) = bc(1)*0.5*dx(1)/k(1) + w(1)
495 IF (bctype(2)) w(n) = bc(2)*0.5*dx(n)/k(n) + w(n)
496 ! print *, 'w(0)=', w(0), 'w(n)=', w(n)
497 ! print *, 'k=', k, 'dx=', dx
498 a = k/dx
499 DO i = 1, n - 1
500 w(i) = (t(i + 1)*a(i + 1) + t(i)*a(i))/(a(i) + a(i + 1))
501 END DO
502 !!FO!!
503 ! PRINT *, 'w: ', w
504 ! PRINT *, 'rhocp: ', rhocp
505 ! PRINT *, 'dx: ', dx
506 ! PRINT *, 'a: ', a
507 DO i = 1, n
508 ! PRINT *, 'i: ', i
509 ! PRINT *, 'i: ', (dt/rhocp(i))
510 ! PRINT *, 'i: ', (w(i - 1) - 2*T(i) + w(i))
511 ! PRINT *, 'i: ', 2*a(i)/dx(i)
512 t1(i) = &
513 (dt/rhocp(i)) &
514 *(w(i - 1) - 2*t(i) + w(i)) &
515 *a(i)/dx(i) &
516 + t(i)
517 END DO
518 !!FO!!
519 ! PRINT *, 'T1: ', T1
520 !for storage the internal distribution of heat should not be important
521 !!FO!! k*d(dT/dx)/dx = rhoCp*(dT/dt) => rhoCp*(dT/dt)*dx = dQs => dQs = k*d(dT/dx)
522 qs = (w(0) - t(1))*2*a(1) + (w(n) - t(n))*2*a(n)
523 ! Qs=sum((T1-T)*rhocp*dx)/dt!
524 t = t1
Referenced by estm_module::estm().
Here is the caller graph for this function:
◆ heatcond1d_ext()
| recursive subroutine heatflux::heatcond1d_ext | ( | real(kind(1d0)), dimension(:), intent(inout) | t, |
| real(kind(1d0)), intent(out) | qs, | ||
| real(kind(1d0)), intent(out) | tsfc, | ||
| real(kind(1d0)), dimension(:), intent(in) | dx, | ||
| real(kind(1d0)), intent(in) | dt, | ||
| real(kind(1d0)), dimension(:), intent(in) | k, | ||
| real(kind(1d0)), dimension(:), intent(in) | rhocp, | ||
| real(kind(1d0)), dimension(2), intent(in) | bc, | ||
| logical, dimension(2), intent(in) | bctype, | ||
| logical, intent(in) | debug ) |
Definition at line 529 of file suews_phys_estm.f95.
530 REAL(KIND(1D0)), INTENT(inout) :: T(:)
531 REAL(KIND(1D0)), INTENT(in) :: dx(:), dt, k(:), rhocp(:), bc(2)
532 REAL(KIND(1D0)), INTENT(out) :: Qs, Tsfc
533 LOGICAL, INTENT(in) :: bctype(2) ! if true, use surrogate flux as boundary condition
534 LOGICAL, INTENT(in) :: debug
535 INTEGER :: i, n !,j !!!!!FO!!!!!
536 REAL(KIND(1D0)), ALLOCATABLE :: w(:), a(:), T1(:), cfl(:)
537 REAL(KIND(1D0)) :: cfl_max
538 REAL(KIND(1D0)), ALLOCATABLE :: T_in(:), T_out(:)
539 REAL(KIND(1D0)) :: dt_x ! for recursion
540 INTEGER :: n_div ! for recursion
541 n = SIZE(t)
542 ALLOCATE (w(0:n), a(n), t1(n), cfl(n), t_in(n), t_out(n))
543
544 ! save initial temperatures
545 t_in = t
546 ! w = interface temperature
547 w(1:n) = t
548 w(0) = bc(1); w(n) = bc(2)
549 !convert from flux to equivalent temperature, not exact
550 ! F = k dT/dX => dx*F/k + T(i+1) = T(i)
551 IF (bctype(1)) w(0) = bc(1)*0.5*dx(1)/k(1) + w(1)
552 IF (bctype(2)) w(n) = bc(2)*0.5*dx(n)/k(n) + w(n)
553 ! print *, 'bc(1)=', bc(1), 'bc(2)=',bc(2)
554 ! print *, 'w(0)=', w(0), 'w(n)=', w(n)
555 ! print *, 'k=', k, 'dx=', dx
556 a = k/dx
557 DO i = 1, n - 1
558 w(i) = (t(i + 1)*a(i + 1) + t(i)*a(i))/(a(i) + a(i + 1))
559 END DO
560 !!FO!!
561 ! PRINT *, 'w: ', w
562 ! PRINT *, 'rhocp: ', rhocp
563 ! PRINT *, 'dt: ', dt
564 ! PRINT *, 'dx: ', dx
565 ! PRINT *, 'a: ', a
566 DO i = 1, n
567 ! PRINT *, 'i: ', i
568 ! PRINT *, 'i: ', (dt/rhocp(i))
569 ! PRINT *, 'i: ', (w(i - 1) - 2*T(i) + w(i))
570 ! PRINT *, 'i: ', 2*a(i)/dx(i)
571 t1(i) = &
572 (dt/rhocp(i)) &
573 *(w(i - 1) - 2*t(i) + w(i)) &
574 *a(i)/dx(i) &
575 + t(i)
576 END DO
577 !!FO!!
578 ! PRINT *, 'T1: ', T1
579 !for storage the internal distribution of heat should not be important
580 !!FO!! k*d(dT/dx)/dx = rhoCp*(dT/dt) => rhoCp*(dT/dt)*dx = dQs => dQs = k*d(dT/dx)
581 ! Qs = (w(0) - T(1))*2*a(1) + (w(n) - T(n))*2*a(n)
582 ! Qs=sum((T1-T)*rhocp*dx)/dt!
583 t = t1
584 tsfc = w(0)
585 ! save output temperatures
586 t_out = t1
587 ! new way for calcualating heat storage
588 qs = sum( &
589 (([bc(1), t_out(1:n - 1)] + t_out)/2. & ! updated temperature
590 -([bc(1), t_in(1:n - 1)] + t_in)/2) & ! initial temperature
591 *rhocp*dx/dt)
592
593 ! fix convergece issue with recursion
594 cfl = abs((w(0:n - 1) - w(1:n))*k/rhocp/(dx**2)*dt)
595 cfl_max = maxval(cfl)
596 ! if (debug) then
597 ! print *, 'T: ', T
598 ! print *, 'rhocp: ', rhocp
599 ! print *, 'w: ', w
600 ! print *, 'a: ', a
601 ! print *, 'cfl: ', cfl
602 ! print *, 'cfl_max: ', cfl_max, MAXLOC(cfl)
603 ! endif
604 IF (cfl_max > .005 .AND. dt > 1) THEN
605 ! re-initialise temperatures at inner nodes
606 t = t_in
607 n_div = 2
608 dt_x = dt/n_div
609 IF (debug) THEN
610 print *, 'entering recursion ', cfl_max, dt_x
611 END IF
612 DO i = 1, n_div
613 CALL heatcond1d_ext(t, qs, tsfc, dx, dt_x, k, rhocp, bc, bctype, debug)
614 END DO
615 END IF
616 ! print *, 'cfl_max: ', cfl_max, MAXLOC(cfl)
References heatcond1d_ext().
Referenced by heatcond1d_ext().
Here is the call graph for this function:
Here is the caller graph for this function:
◆ heatcond1d_vstep()
| subroutine heatflux::heatcond1d_vstep | ( | real(kind(1d0)), dimension(:), intent(inout) | t, |
| real(kind(1d0)), intent(out) | qs, | ||
| real(kind(1d0)), intent(out) | tsfc, | ||
| real(kind(1d0)), dimension(:), intent(in) | dx, | ||
| real(kind(1d0)), intent(in) | dt, | ||
| real(kind(1d0)), dimension(:), intent(in) | k, | ||
| real(kind(1d0)), dimension(:), intent(in) | rhocp, | ||
| real(kind(1d0)), dimension(2), intent(in) | bc, | ||
| logical, dimension(2), intent(in) | bctype, | ||
| logical, intent(in) | debug ) |
Definition at line 4 of file suews_phys_ehc.f95.
5 REAL(KIND(1D0)), INTENT(inout) :: T(:)
6 REAL(KIND(1D0)), INTENT(in) :: dx(:), dt, k(:), rhocp(:), bc(2)
7 REAL(KIND(1D0)), INTENT(out) :: Qs, Tsfc
8 LOGICAL, INTENT(in) :: bctype(2) ! if true, use surrogate flux as boundary condition
9 LOGICAL, INTENT(in) :: debug
10 INTEGER :: i, n
11 REAL(KIND(1D0)), ALLOCATABLE :: w(:), a(:), T1(:), cfl(:)
12
13 REAL(KIND(1D0)) :: cfl_max
14 REAL(KIND(1D0)) :: dt_remain
15 REAL(KIND(1D0)) :: dt_step
16 REAL(KIND(1D0)) :: dt_step_cfl
17
18 REAL(KIND(1D0)), ALLOCATABLE :: T_in(:), T_out(:)
19 REAL(KIND(1D0)) :: dt_x ! for recursion
20 INTEGER :: n_div ! for recursion
21 n = SIZE(t)
22 ALLOCATE (w(0:n), a(n), t1(n), cfl(n), t_in(n), t_out(n))
23
24 ! save initial temperatures
25 t_in = t
26 ! w = interface temperature
27 w(1:n) = t
28 w(0) = bc(1); w(n) = bc(2)
29 !convert from flux to equivalent temperature, not exact
30 ! F = k dT/dX => dx*F/k + T(i+1) = T(i)
31 IF (bctype(1)) w(0) = bc(1)*0.5*dx(1)/k(1) + w(1)
32 IF (bctype(2)) w(n) = bc(2)*0.5*dx(n)/k(n) + w(n)
33 ! print *, 'bc(1)=', bc(1), 'bc(2)=',bc(2)
34 ! print *, 'w(0)=', w(0), 'w(n)=', w(n)
35 ! print *, 'k=', k, 'dx=', dx
36 a = k/dx
37 DO i = 1, n - 1
38 w(i) = (t(i + 1)*a(i + 1) + t(i)*a(i))/(a(i) + a(i + 1))
39 END DO
40
41 ! fix convergece issue with recursion
42 dt_remain = dt
43 dt_step_cfl = 0.05*minval(dx**2/(k/rhocp))
44 DO WHILE (dt_remain > 1e-10)
45 dt_step = min(dt_step_cfl, dt_remain)
46 !PRINT *, 'dt_remain: ', dt_remain
47 !PRINT *, 'dt_step_cfl: ', dt_step_cfl
48 !PRINT *, '***** dt_step: ', dt_step
49 !!FO!!
50 ! PRINT *, 'w: ', w
51 ! PRINT *, 'rhocp: ', rhocp
52 ! PRINT *, 'dt: ', dt
53 ! PRINT *, 'dx: ', dx
54 ! PRINT *, 'a: ', a
55 DO i = 1, n
56 ! PRINT *, 'i: ', i
57 ! PRINT *, 'i: ', (dt/rhocp(i))
58 ! PRINT *, 'i: ', (w(i - 1) - 2*T(i) + w(i))
59 ! PRINT *, 'i: ', 2*a(i)/dx(i)
60 t1(i) = &
61 (dt_step/rhocp(i)) &
62 *(w(i - 1) - 2*t(i) + w(i)) &
63 *a(i)/dx(i) &
64 + t(i)
65 END DO
66 t = t1
67 DO i = 1, n - 1
68 w(i) = (t(i + 1)*a(i + 1) + t(i)*a(i))/(a(i) + a(i + 1))
69 END DO
70 dt_remain = dt_remain - dt_step
71 END DO
72 !!FO!!
73 ! PRINT *, 'T1: ', T1
74 !for storage the internal distribution of heat should not be important
75 !!FO!! k*d(dT/dx)/dx = rhoCp*(dT/dt) => rhoCp*(dT/dt)*dx = dQs => dQs = k*d(dT/dx)
76 ! Qs = (w(0) - T(1))*2*a(1) + (w(n) - T(n))*2*a(n)
77 ! Qs=sum((T1-T)*rhocp*dx)/dt!
78
79 tsfc = w(0)
80 ! save output temperatures
81 t_out = t1
82 ! new way for calcualating heat storage
83 qs = sum( &
84 (([bc(1), t_out(1:n - 1)] + t_out)/2. & ! updated temperature
85 -([bc(1), t_in(1:n - 1)] + t_in)/2) & ! initial temperature
86 *rhocp*dx/dt)
Referenced by ehc_module::ehc().
Here is the caller graph for this function:
