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 solves the vertical transport equation.
6 : !! <cvert> is temporary storage for concentrations (particles,
7 : !! gas, potential temperature) being transported.
8 : !! New values of <cvert> are calculated.
9 : !!
10 : !! @author Eric Jensen
11 : !! @version Dec-1996
12 231137280 : subroutine versol(carma, cstate, cvert, itbnd, ibbnd, ftop, fbot, cvert_tbnd, cvert_bbnd, &
13 231137280 : vertadvu, vertadvd, vertdifu, vertdifd, rc)
14 :
15 : ! types
16 : use carma_precision_mod
17 : use carma_enums_mod
18 : use carma_constants_mod
19 : use carma_types_mod
20 : use carmastate_mod
21 : use carma_mod
22 :
23 : implicit none
24 :
25 : type(carma_type), intent(in) :: carma !! the carma object
26 : type(carmastate_type), intent(inout) :: cstate !! the carma state object
27 : real(kind=f), intent(inout) :: cvert(NZ) !! quantity being transported
28 : integer, intent(in) :: itbnd !! top boundary condition
29 : integer, intent(in) :: ibbnd !! bottom boundary condition
30 : real(kind=f), intent(in) :: ftop !! flux at top boundary
31 : real(kind=f), intent(in) :: fbot !! flux at bottom boundary
32 : real(kind=f), intent(in) :: cvert_tbnd !! quantity at top boundary
33 : real(kind=f), intent(in) :: cvert_bbnd !! quantity at bottom boundary
34 : real(kind=f), intent(in) :: vertadvu(NZP1) !! upward vertical transport rate into level k from level k-1 [cm/s]
35 : real(kind=f), intent(in) :: vertadvd(NZP1) !! downward vertical transport rate into level k from level k-1 [cm/s]
36 : real(kind=f), intent(in) :: vertdifu(NZP1) !! upward vertical diffusion rate into level k from level k-1 [cm/s]
37 : real(kind=f), intent(in) :: vertdifd(NZP1) !! downward vertical diffusion rate into level k from level k-1 [cm/s]
38 : integer, intent(inout) :: rc !! return code, negative indicates failure
39 :
40 : ! Declare local variables
41 : !
42 : integer :: k
43 462274560 : real(kind=f) :: al(NZ)
44 462274560 : real(kind=f) :: bl(NZ)
45 462274560 : real(kind=f) :: dl(NZ)
46 462274560 : real(kind=f) :: el(NZ)
47 462274560 : real(kind=f) :: fl(NZ)
48 462274560 : real(kind=f) :: ul(NZ)
49 462274560 : real(kind=f) :: ctempl(NZ)
50 462274560 : real(kind=f) :: ctempu(NZ)
51 231137280 : real(kind=f) :: divcor(NZ)
52 : real(kind=f) :: uc
53 : real(kind=f) :: cour
54 : real(kind=f) :: denom
55 :
56 : ! Divergence adjustments are not being generated.
57 7627530240 : divcor(:) = 0._f
58 :
59 : ! Determine whether transport should be solved explicitly (uc=0)
60 : ! or implicitly (uc=1).
61 231137280 : uc = 0._f
62 7627530240 : do k = 1,NZ
63 0 : cour = dz(k)/dtime - &
64 7396392960 : ( vertdifu(k+1) + vertdifd(k) + vertadvu(k+1) + vertadvd(k) )
65 :
66 7627530240 : if( cour .lt. 0._f .and. uc .ne. 1._f )then
67 11166406 : uc = 1.0_f
68 :
69 : ! NOTE: This can happen a lot and clutters up the log. Should we print it out or not?
70 : ! write(LUNOPRT,'(a,i3,7(1x,1pe8.1))') &
71 : ! 'in versol: k dz/dt vdifd vdifu vadvd vadvu cour uc = ', &
72 : ! k, dz(k)/dtime, vertdifd(k), vertdifu(k+1), &
73 : ! vertadvd(k), vertadvu(k+1), cour, uc
74 : endif
75 : enddo
76 :
77 : ! Store concentrations in local variables (shifted up and down
78 : ! a vertical level).
79 7396392960 : do k = 2,NZ
80 7165255680 : ctempl(k) = cvert(k-1)
81 7396392960 : ctempu(k-1) = cvert(k)
82 : enddo
83 :
84 231137280 : if( ibbnd .eq. I_FIXED_CONC ) then
85 231137280 : ctempl(1) = cvert_bbnd
86 : else
87 0 : ctempl(1) = 0._f
88 : endif
89 :
90 231137280 : if( itbnd .eq. I_FIXED_CONC ) then
91 231137280 : ctempu(NZ) = cvert_tbnd
92 : else
93 0 : ctempu(NZ) = 0._f
94 : endif
95 :
96 : ! Calculate coefficients of the transport equation:
97 : ! al(k)c(k+1) + bl(k)c(k) + ul(k)c(k-1) = dl(k)
98 :
99 7627530240 : do k = 1,NZ
100 7396392960 : al(k) = uc * ( vertdifd(k+1) + vertadvd(k+1) )
101 7396392960 : bl(k) = -( uc*(vertdifd(k)+vertdifu(k+1)+ &
102 : vertadvd(k)+vertadvu(k+1)) &
103 7396392960 : + dz(k)/dtime )
104 7396392960 : ul(k) = uc * ( vertdifu(k) + vertadvu(k) )
105 : dl(k) = cvert(k) * &
106 : ( (1._f - uc)*(vertdifd(k)+vertdifu(k+1)+ &
107 : vertadvd(k)+vertadvu(k+1)) &
108 0 : - dz(k)/dtime ) - &
109 : (1._f - uc) * ( (vertdifu(k)+vertadvu(k))*ctempl(k) + &
110 : (vertdifd(k+1)+vertadvd(k+1))*ctempu(k) ) - &
111 7627530240 : divcor(k) * dz(k)
112 : enddo
113 :
114 : ! Boundary fluxes: <ftop> is the downward flux across the
115 : ! upper boundary; <fbot> is the upward flux across the
116 : ! lower boundary.
117 231137280 : if(( igridv .eq. I_SIG ) .or. ( igridv .eq. I_HYBRID )) then
118 231137280 : if( itbnd .eq. I_FLUX_SPEC ) dl(1) = dl(1) - ftop
119 231137280 : if( ibbnd .eq. I_FLUX_SPEC ) dl(NZ) = dl(NZ) - fbot
120 : else
121 0 : if( itbnd .eq. I_FLUX_SPEC ) dl(NZ) = dl(NZ) - ftop
122 0 : if( ibbnd .eq. I_FLUX_SPEC ) dl(1) = dl(1) - fbot
123 : endif
124 :
125 : ! Calculate recursion relations.
126 231137280 : el(1) = dl(1)/bl(1)
127 231137280 : fl(1) = al(1)/bl(1)
128 7396392960 : do k = 2,NZ
129 7165255680 : denom = bl(k) - ul(k) * fl(k-1)
130 7165255680 : el(k) = ( dl(k) - ul(k)*el(k-1) ) / denom
131 7396392960 : fl(k) = al(k) / denom
132 : enddo
133 :
134 : ! Calculate new concentrations.
135 :
136 231137280 : cvert(NZ) = el(NZ)
137 7396392960 : do k = NZ-1,1,-1
138 7396392960 : cvert(k) = el(k) - fl(k)*cvert(k+1)
139 : enddo
140 :
141 : ! Return to caller with new concentrations.
142 231137280 : return
143 231137280 : end
|
