### File: C:\NOAA\NEMS_11731\src\chem\gocart\src\Components\GOCART_GridComp\CARMA_GridComp\versol.F90

```1     ! Colarco, May 14, 2007
2     ! sed command to replace F77 comments with F90: sed 's/^c/\!/g'
3     ! F90-version of original CARMA versol.f routine (see comments below from
5     !
6     ! Notes:
7     !  1) I don't allow divcor for now.
8     !
9        subroutine versol ( ix, iy, km, ibbnd, itbnd, fbot, ftop, &
11                            cvert_bbnd, cvert_tbnd, cvert, carma, rc )
12     !
13     ! Define variables and usages
14        use carma_types_mod
15
16        implicit none
17
18     !  Input/Output
19        integer, intent(out)             :: rc         ! Error return code:
20                                                       !  0 - all is well
21                                                       !  1 -
22     !  Input/Output
23        integer :: ix, iy, km, itbnd, ibbnd
24        real(kind=f) :: fbot, ftop
26                                         vertdifu, vertdifd
27        real(kind=f), dimension(km)   :: cvert
28        real(kind=f) :: cvert_bbnd, cvert_tbnd
29
30     !  Local
31        integer :: k
32        real(kind=f) :: uc, cour, denom
33        real(kind=f), dimension(km)   :: al, bl, dl, ul, el, fl, &
34                                         ctempl, ctempu, divcor
35
36     #include "carma_globaer.h"
37
38     #ifdef DEBUG
39        write(*,*) '+ versol'
40     #endif
41
42        divcor(:) = 0._f
43
44     !
45     !
46     !  @(#) versol.f  Jensen  Dec-1996
47     !  This routine solves the vertical transport equation.
48     !  <cvert> is temporary storage for concentrations (particles,
49     !  gas, potential temperature) being transported.
50     !  New values of <cvert> are calculated.
51     !
52     !  Modified  Sep-1997  (McKie)
53     !  Remove <ixy> from arg list
54     !  <ixy> now available as a global var in common block.
55     !
56     !  Argument list input:
57     !    None.
58     !
59     !  Argument list output:
60     !    None.
61     !
62     !
63     !  Include global constants and variables.
64     !
65     !      include 'globaer.h'
66     !
67     !
68     !  Declare local variables
69     !
70     !      dimension al(NZ), bl(NZ), dl(NZ), ul(NZ),
71     !     \$          el(NZ), fl(NZ), ctempl(NZ), ctempu(NZ)
72     !
73     !
74     !  Announce entry to this routine.
75     !
76     !      if( DEBUG ) write(LUNOPRT,'(/,a)') 'Enter versol'
77     !
78     !
79     !  Determine whether transport should be solved explicitly (uc=0)
80     !  or implicitly (uc=1).
81     !
82           uc = 0._f
83           do k = 1,NZ
84             cour = dz(ix,iy,k)/dtime - &
86             if( cour .lt. 0. .and. uc .ne. 1. )then
87               uc = 1.0_f
88     !          write(LUNOPRT,'(a,i3,7(1x,1pe8.1))')
89     !     \$      'in versol: k dz/dt vdifd vdifu vadvd vadvu cour uc = ',
90     !     \$       k, dz(ix,iy,k)/dtime, vertdifd(k), vertdifu(k+1),
92             endif
93           enddo
94     !
95     !
96     !  Store concentrations in local variables (shifted up and down
97     !  a vertical level).
98     !
99           do k = 2,NZ
100             ctempl(k) = cvert(k-1)
101             ctempu(k-1) = cvert(k)
102           enddo
103
104           if( ibbnd .eq. I_FIXED_CONC ) then
105             ctempl(1) = cvert_bbnd
106           else
107             ctempl(1) = 0._f
108           endif
109
110           if( itbnd .eq. I_FIXED_CONC ) then
111             ctempu(NZ) = cvert_tbnd
112           else
113            ctempu(NZ) = 0._f
114           endif
115     !
116     !
117     !  Calculate coefficients of the transport equation:
118     !    al(k)c(k+1) + bl(k)c(k) + ul(k)c(k-1) = dl(k)
119     !  This part could use a notation.  The equation solved is really
120     !  a very simply formulated solution to the flux form of the
122     !    du/dt = d(vu)/dx
123     !  This is solved explicitly if the cfl condition is met:
124     !
125     !     n+1    n            n      n                n    n
126     !    u    - u  = v      (u    - u  )  +  v      (u  - u   )
127     !     k      k    k+1/2   k+1    k        k-1/2   k    k-1
128     !    ---------   ------------------------------------------
129     !       dt                           dx
130     !
131     !  or solved implicitly if it is not:
132     !
133     !     n+1    n            n+1    n+1              n+1  n+1
134     !    u    - u  = v      (u    - u  )  +  v      (u  - u   )
135     !     k      k    k+1/2   k+1    k        k-1/2   k    k-1
136     !    ---------   ------------------------------------------
137     !       dt                           dx
138     !
139     !  Where v(k+1/2) = vertadvd(k+1)+vertdifd(k+1)
140     !    and v(k-1/2) = vertadvu(k)+vertdifu(k)
141     !  The n+1 time-level terms are the left side of the equation and the
142     !  n time-level terms are the right side (we know the right side exactly).
143     !  The switch u is used to flip between implicit and explicit solutions.
144     !  The resulting set of coefficients are put into a tridiagonal solver
145     !  (below) with possibly a flux boundary condition implied.
146     !
147           do k = 1,NZ
148             al(k) = uc * ( vertdifd(k+1) + vertadvd(k+1) )
149             bl(k) = -( uc*(vertdifd(k)+vertdifu(k+1)+ &
151                        + dz(ix,iy,k)/dtime )
152             ul(k) = uc * ( vertdifu(k) + vertadvu(k) )
153             dl(k) = cvert(k) * &
154                ( (1.-uc)*(vertdifd(k)+vertdifu(k+1)+ &
156                - dz(ix,iy,k)/dtime ) - &
157                (1.-uc) * ( (vertdifu(k)+vertadvu(k))*ctempl(k) + &
159                divcor(k) * dz(ix,iy,k)
160           enddo
161     !
162     !
163     !  Boundary fluxes: <ftop> is the downward flux across the
164     !  upper boundary; <fbot> is the upward flux across the
165     !  lower boundary.  Top and bottom in the coordinate sense.
166     !
167            if( itbnd .eq. I_FLUX_SPEC ) dl(NZ) = dl(NZ) - ftop
168            if( ibbnd .eq. I_FLUX_SPEC ) dl(1) = dl(1) - fbot
169     !
170     !
171     !  Calculate recursion relations.
172     !
173           el(1) = dl(1)/bl(1)
174           fl(1) = al(1)/bl(1)
175           do k = 2,NZ
176             denom = bl(k) - ul(k) * fl(k-1)
177             el(k) = ( dl(k) - ul(k)*el(k-1) ) / denom
178             fl(k) = al(k) / denom
179           enddo
180     !
181     !
182     !  Calculate new concentrations.
183     !
184           cvert(NZ) = el(NZ)
185           do k = NZ-1,1,-1
186             cvert(k) = el(k) - fl(k)*cvert(k+1)
187           enddo
188     !
189     !