Line data Source code
1 : ! Include shortname defintions, so that the F77 code does not have to be modified to
2 : ! reference the CARMA structure.
3 : #include "carma_globaer.h"
4 :
5 : !! This routine calculates new potential temperature concentration
6 : !! (and updates temperature) due to microphysical and radiative forcings.
7 : !! The equation solved (the first law of thermodynamics in flux form) is
8 : !!
9 : !! d(rhostar theta) rhostar theta d(qv) 1 dF
10 : !! --------------- = - ------------- * ( L ----- + --- -- )
11 : !! dt Cp T dt rho dz
12 : !!
13 : !! where
14 : !! rhostar = scaled air density
15 : !! theta = potential temperature
16 : !! t = time
17 : !! Cp = specific heat (at constant pressure) of air
18 : !! T = air temperature
19 : !! qv = water vapor mixing ratio
20 : !! L = latent heat
21 : !! F = net radiative flux
22 : !! z = unscaled altitude
23 : !!
24 : !! @author Andy Ackerman
25 : !! @version Oct-1997
26 0 : subroutine tsolve(carma, cstate, iz, scale_threshold, rc)
27 :
28 : ! types
29 : use carma_precision_mod
30 : use carma_enums_mod
31 : use carma_constants_mod
32 : use carma_types_mod
33 : use carmastate_mod
34 : use carma_mod
35 :
36 : implicit none
37 :
38 : type(carma_type), intent(in) :: carma !! the carma object
39 : type(carmastate_type), intent(inout) :: cstate !! the carma state object
40 : integer, intent(in) :: iz !! z index
41 : real(kind=f) :: scale_threshold !! Scaling factor for convergence thresholds
42 : integer, intent(inout) :: rc !! return code, negative indicates failure
43 :
44 : 1 format(/,'tsolve::ERROR - negative temperature for : iz=',i4,',xc=',&
45 : f7.2,',yc=',f7.2,',T=',e10.3,',dT=',e10.3,',t_old=',e10.3,',d_gc=',e10.3,',dT_adv=',e10.3)
46 : 2 format(/,'tsolve::ERROR - temperature change to large for : iz=',i4,',xc=',&
47 : f7.2,',yc=',f7.2,',T=',e10.3,',dT_rlh=',e10.3,',dT_pth=',e10.3,',t_old=',e10.3,',d_gc=',e10.3,',dT_adv=',e10.3)
48 : 3 format(/,'tsolve::ERROR - temperature change to large for : iz=',i4,',xc=',&
49 : f7.2,',yc=',f7.2,',T=',e10.3,',dT_rlh=',e10.3,',dT_pth=',e10.3,',t_old=',e10.3)
50 : 4 format(/,'tsolve::ERROR - negative temperature for : iz=',i4,',xc=',&
51 : f7.2,',yc=',f7.2,',T=',e10.3,',dT=',e10.3,',t_old=',e10.3)
52 :
53 : real(kind=f) :: dt ! delta temperature
54 : real(kind=f) :: threshold ! convergence threshold
55 :
56 :
57 : ! Solve for the new <t> due to latent heat exchange and radiative heating.
58 : ! Latent and radiative heating rates are in units of [deg_K/s].
59 : !
60 : ! NOTE: In the embedded model rhoa and p are handled by the parent model and
61 : ! won't change during one time step.
62 : !
63 : ! NOTE: Radiative heating by the particles is handled by the parent model, so
64 : ! that term does not need to be added here.
65 0 : dt = dtime * rlprod
66 0 : rlheat(iz) = rlheat(iz) + rlprod * dtime
67 :
68 : ! With particle heating, you must also include the impact of heat
69 : ! conduction from the particle
70 : !
71 : ! NOTE: We are ignoring the energy to heat the particle, since we
72 : ! are not tracking the particle temperature. Thus ...
73 0 : if (do_pheatatm) then
74 0 : dt = dt + dtime * phprod
75 0 : partheat(iz) = partheat(iz) + phprod * dtime
76 : end if
77 :
78 0 : t(iz) = t(iz) + dt
79 :
80 :
81 : ! Don't let the temperature go negative.
82 0 : if (t(iz) < 0._f) then
83 0 : if (do_substep) then
84 0 : if (nretries == maxretries) then
85 0 : if (do_print) write(LUNOPRT,1) iz, xc, yc, t(iz), dt, told(iz), d_gc(iz, 1), d_t(iz)
86 : end if
87 : else
88 0 : if (do_print) write(LUNOPRT,4) iz, xc, yc, t(iz), dt, told(iz)
89 : end if
90 :
91 0 : rc = RC_WARNING_RETRY
92 : end if
93 :
94 : ! Don't let the temperature change by more than the threshold in any given substep,
95 : ! to prevent overshooting that doesn't result in negative gas concentrations, but
96 : ! does result in excessive temperature swings.
97 0 : threshold = dt_threshold / scale_threshold
98 :
99 0 : if (threshold /= 0._f) then
100 0 : if (abs(abs(dt)) > threshold) then
101 0 : if (do_substep) then
102 0 : if (nretries == maxretries) then
103 0 : if (do_print) write(LUNOPRT,2) iz, xc, yc, t(iz), rlprod*dtime, dtime*partheat(iz), told(iz), d_gc(iz, 1), d_t(iz)
104 : end if
105 : else
106 0 : if (do_print) write(LUNOPRT,3) iz, xc, yc, t(iz), rlprod*dtime, dtime*partheat(iz), told(iz)
107 : end if
108 :
109 0 : rc = RC_WARNING_RETRY
110 : end if
111 : end if
112 :
113 : ! Return to caller with new temperature.
114 0 : return
115 0 : end
|
