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 4 ! original routine header). 5 ! 6 ! Notes: 7 ! 1) I don't allow divcor for now. 8 ! 9 subroutine versol ( ix, iy, km, ibbnd, itbnd, fbot, ftop, & 10 vertadvu, vertadvd, vertdifu, vertdifd, & 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 25 real(kind=f), dimension(km+1) :: vertadvu, vertadvd, & 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 - & 85 ( vertdifu(k+1) + vertdifd(k) + vertadvu(k+1) + vertadvd(k) ) 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), 91 ! $ vertadvd(k), vertadvu(k+1), cour, uc 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 121 ! advection equation: 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)+ & 150 vertadvd(k)+vertadvu(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)+ & 155 vertadvd(k)+vertadvu(k+1)) & 156 - dz(ix,iy,k)/dtime ) - & 157 (1.-uc) * ( (vertdifu(k)+vertadvu(k))*ctempl(k) + & 158 (vertdifd(k+1)+vertadvd(k+1))*ctempu(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 ! 190 ! Return to caller with new concentrations. 191 ! 192 rc = 0 193 194 return 195 end 196