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 the number of sub-timesteps <ntsubsteps>
6 : !! for the current model spatial point.
7 : !!
8 : !! @author Eric Jensen
9 : !! @version Apr-2000
10 73543680 : subroutine nsubsteps(carma, cstate, iz, dtime_save, ntsubsteps, rc)
11 :
12 : ! types
13 : use carma_precision_mod
14 : use carma_enums_mod
15 : use carma_constants_mod
16 : use carma_types_mod
17 : use carmastate_mod
18 : use carma_mod
19 :
20 : implicit none
21 :
22 : type(carma_type), intent(in) :: carma !! the carma object
23 : type(carmastate_type), intent(inout) :: cstate !! the carma state object
24 : integer, intent(in) :: iz !! z index
25 : real(kind=f), intent(in) :: dtime_save !! original (not substepped) dtime
26 : integer, intent(inout) :: ntsubsteps !! suggested number of substeps
27 : integer, intent(inout) :: rc !! return code, negative indicates failure
28 :
29 : ! Local declarations
30 : integer :: ig ! group index
31 : integer :: igas ! gas index
32 : integer :: ibin ! bin index
33 : integer :: iepart
34 : integer :: inuc
35 : integer :: ienucto
36 147087360 : integer :: ibin_small(NGROUP)
37 : real(kind=f) :: g0
38 : real(kind=f) :: g1
39 : real(kind=f) :: dmdt
40 : real(kind=f) :: dt_adv
41 : real(kind=f) :: ss
42 : real(kind=f) :: ssold
43 : real(kind=f) :: pvap
44 : real(kind=f) :: vf_max
45 :
46 :
47 : ! If substepping is disabled, then use one substep
48 73543680 : if (.not. do_substep) then
49 0 : ntsubsteps = 1
50 : else
51 : ! Set default values
52 73543680 : ntsubsteps = minsubsteps
53 :
54 : ! Find the bin number of the smallest particle bin that
55 : ! contains a significant number of particles.
56 : ! Also check for significant activation of water droplets.
57 :
58 73543680 : if( ntsubsteps .lt. maxsubsteps )then
59 :
60 220631040 : do ig = 1, NGROUP
61 :
62 220631040 : if( pconmax(iz,ig) .gt. FEW_PC) then
63 :
64 46313084 : ibin_small(ig) = NBIN
65 :
66 : ! element of particle number concentration
67 46313084 : iepart = ienconc(ig)
68 :
69 46313084 : if( itype(iepart) .eq. I_INVOLATILE ) then
70 :
71 : ! condensing gas
72 0 : igas = inucgas(ig)
73 :
74 0 : if (igas /= 0) then
75 :
76 0 : ss = max( supsatl(iz,igas), supsatlold(iz,igas) )
77 :
78 0 : do inuc = 1,nnuc2elem(iepart)
79 0 : ienucto = inuc2elem(inuc,iepart)
80 :
81 0 : if( inucproc(iepart,ienucto) .eq. I_DROPACT ) then
82 0 : do ibin = 1, NBIN
83 0 : if( pc(iz,ibin,iepart) / zmet(iz) .gt. conmax * pconmax(iz,ig) .and. &
84 0 : ss .gt. scrit(iz,ibin,ig) )then
85 0 : ntsubsteps = maxsubsteps
86 : endif
87 : enddo
88 : endif
89 : enddo
90 : endif
91 :
92 46313084 : elseif( itype(iepart) .eq. I_VOLATILE ) then
93 :
94 926261680 : do ibin = NBIN-1, 1, -1
95 926261680 : if( pc(iz,ibin,iepart) / zmet(iz) .gt. conmax * pconmax(iz,ig) )then
96 139640441 : ibin_small(ig) = ibin
97 : endif
98 : enddo
99 :
100 : endif
101 : endif
102 : enddo
103 : endif
104 :
105 : ! Calculate the growth rate of a particle with the mode radius for
106 : ! each volatile group. The maximum time-step to use is then the
107 : ! mass growth rate divided by the mass bin width / 2.
108 73543680 : if( ntsubsteps .lt. maxsubsteps )then
109 :
110 73543680 : dt_adv = dtime_save
111 220631040 : do ig = 1, NGROUP
112 :
113 : ! element of particle number concentration
114 147087360 : iepart = ienconc(ig)
115 :
116 : ! condensing gas
117 147087360 : igas = igrowgas(iepart)
118 :
119 220631040 : if (igas /= 0) then
120 :
121 147087360 : if( pconmax(iz,ig) .gt. FEW_PC ) then
122 :
123 46313084 : if( itype(iepart) .eq. I_VOLATILE ) then
124 :
125 46313084 : if( is_grp_ice(ig) )then
126 0 : ss = supsati(iz,igas)
127 0 : pvap = pvapi(iz,igas)
128 : else
129 46313084 : ss = supsatl(iz,igas)
130 46313084 : pvap = pvapl(iz,igas)
131 : endif
132 :
133 46313084 : g0 = gro(iz,ibin_small(ig),ig)
134 46313084 : g1 = gro1(iz,ibin_small(ig),ig)
135 46313084 : dmdt = abs( pvap * ss * g0 / ( 1._f + g0*g1*pvap ) )
136 :
137 46313084 : if (dmdt /= 0._f) then
138 46313084 : dt_adv = min( dt_adv, dm(ibin_small(ig),ig)/dmdt )
139 : end if
140 : endif
141 : endif
142 : endif
143 : enddo
144 :
145 73543680 : ntsubsteps = nint(min(real(maxsubsteps, kind=f), real(dtime_save, kind=f) / dt_adv))
146 73543680 : ntsubsteps = max( minsubsteps, ntsubsteps )
147 : endif
148 :
149 : ! If the ice supersaturation is large enough for homogeneous freezing
150 : ! of sulfate aerosols, then use maximum number of substeps
151 73543680 : if( ntsubsteps .lt. (maxsubsteps) )then
152 220548696 : do ig = 1, NGROUP
153 :
154 : ! element of particle number concentration
155 147032464 : iepart = ienconc(ig)
156 :
157 : ! condensing gas
158 147032464 : igas = inucgas(ig)
159 :
160 220548696 : if (igas /= 0) then
161 :
162 147032464 : do inuc = 1,nnuc2elem(iepart)
163 73516232 : ienucto = inuc2elem(inuc,iepart)
164 :
165 147032464 : if (iand(inucproc(iepart,ienucto), I_AERFREEZE) .ne. 0) then
166 0 : if( (supsati(iz,igas) .gt. 0.4_f) .and. (t(iz) .lt. 233.16_f) ) then
167 0 : ntsubsteps = maxsubsteps
168 : endif
169 : endif
170 : enddo
171 : endif
172 : enddo
173 : endif
174 : endif
175 :
176 : ! Return to caller with number of sub-timesteps evaluated.
177 73543680 : return
178 73543680 : end
|
