| 1 | !********************************************************************** |
|---|
| 2 | ! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 * |
|---|
| 3 | ! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, * |
|---|
| 4 | ! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann * |
|---|
| 5 | ! * |
|---|
| 6 | ! This file is part of FLEXPART. * |
|---|
| 7 | ! * |
|---|
| 8 | ! FLEXPART is free software: you can redistribute it and/or modify * |
|---|
| 9 | ! it under the terms of the GNU General Public License as published by* |
|---|
| 10 | ! the Free Software Foundation, either version 3 of the License, or * |
|---|
| 11 | ! (at your option) any later version. * |
|---|
| 12 | ! * |
|---|
| 13 | ! FLEXPART is distributed in the hope that it will be useful, * |
|---|
| 14 | ! but WITHOUT ANY WARRANTY; without even the implied warranty of * |
|---|
| 15 | ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * |
|---|
| 16 | ! GNU General Public License for more details. * |
|---|
| 17 | ! * |
|---|
| 18 | ! You should have received a copy of the GNU General Public License * |
|---|
| 19 | ! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. * |
|---|
| 20 | !********************************************************************** |
|---|
| 21 | |
|---|
| 22 | subroutine redist (ipart,ktop,ipconv) |
|---|
| 23 | |
|---|
| 24 | !************************************************************************** |
|---|
| 25 | ! Do the redistribution of particles due to convection |
|---|
| 26 | ! This subroutine is called for each particle which is assigned |
|---|
| 27 | ! a new vertical position randomly, based on the convective redistribution |
|---|
| 28 | ! matrix |
|---|
| 29 | !************************************************************************** |
|---|
| 30 | |
|---|
| 31 | ! Petra Seibert, Feb 2001, Apr 2001, May 2001, Jan 2002, Nov 2002 and |
|---|
| 32 | ! Andreas Frank, Nov 2002 |
|---|
| 33 | |
|---|
| 34 | ! Caroline Forster: November 2004 - February 2005 |
|---|
| 35 | |
|---|
| 36 | use par_mod, only: ga, r_air, nuvzmax |
|---|
| 37 | use com_mod, only: ztra1, ldirect, lsynctime, nz, height |
|---|
| 38 | use conv_mod, only: fmassfrac, nconvtop, pconv, phconv, dpr, tconv, qconv, psconv, & |
|---|
| 39 | td2conv, tt2conv, sub, uvzlev |
|---|
| 40 | use random_mod, only: ran3_conv |
|---|
| 41 | |
|---|
| 42 | implicit none |
|---|
| 43 | |
|---|
| 44 | real,parameter :: const=r_air/ga |
|---|
| 45 | integer :: ipart, ktop,ipconv |
|---|
| 46 | integer :: k, kz, levnew, levold |
|---|
| 47 | real :: wsub(nuvzmax) |
|---|
| 48 | real :: totlevmass, wsubpart |
|---|
| 49 | real :: temp_levold,temp_levold1 |
|---|
| 50 | real :: sub_levold,sub_levold1 |
|---|
| 51 | real :: pint, pold, rn, tv, tvold, dlevfrac |
|---|
| 52 | real :: ew,ztold,ffraction |
|---|
| 53 | real :: tv1, tv2, dlogp, dz, dz1, dz2 |
|---|
| 54 | |
|---|
| 55 | |
|---|
| 56 | ! ipart ... number of particle to be treated |
|---|
| 57 | |
|---|
| 58 | ipconv=1 |
|---|
| 59 | |
|---|
| 60 | ! determine height of the eta half-levels (uvzlev) |
|---|
| 61 | ! do that only once for each grid column |
|---|
| 62 | ! i.e. when ktop.eq.1 |
|---|
| 63 | !************************************************************** |
|---|
| 64 | |
|---|
| 65 | if (ktop .le. 1) then |
|---|
| 66 | |
|---|
| 67 | tvold=tt2conv*(1.+0.378*ew(td2conv)/psconv) |
|---|
| 68 | pold=psconv |
|---|
| 69 | uvzlev(1)=0. |
|---|
| 70 | |
|---|
| 71 | pint = phconv(2) |
|---|
| 72 | ! determine next virtual temperatures |
|---|
| 73 | tv1 = tconv(1)*(1.+0.608*qconv(1)) |
|---|
| 74 | tv2 = tconv(2)*(1.+0.608*qconv(2)) |
|---|
| 75 | ! interpolate virtual temperature to half-level |
|---|
| 76 | tv = tv1 + (tv2-tv1)*(pconv(1)-phconv(2))/(pconv(1)-pconv(2)) |
|---|
| 77 | if (abs(tv-tvold).gt.0.2) then |
|---|
| 78 | uvzlev(2) = uvzlev(1) + & |
|---|
| 79 | const*log(pold/pint)* & |
|---|
| 80 | (tv-tvold)/log(tv/tvold) |
|---|
| 81 | else |
|---|
| 82 | uvzlev(2) = uvzlev(1)+ & |
|---|
| 83 | const*log(pold/pint)*tv |
|---|
| 84 | endif |
|---|
| 85 | tvold=tv |
|---|
| 86 | tv1=tv2 |
|---|
| 87 | pold=pint |
|---|
| 88 | |
|---|
| 89 | ! integrate profile (calculation of height agl of eta layers) as required |
|---|
| 90 | do kz = 3, nconvtop+1 |
|---|
| 91 | ! note that variables defined in calcmatrix.f (pconv,tconv,qconv) |
|---|
| 92 | ! start at the first real ECMWF model level whereas kz and |
|---|
| 93 | ! thus uvzlev(kz) starts at the surface. uvzlev is defined at the |
|---|
| 94 | ! half-levels (between the tconv, qconv etc. values !) |
|---|
| 95 | ! Thus, uvzlev(kz) is the lower boundary of the tconv(kz) cell. |
|---|
| 96 | pint = phconv(kz) |
|---|
| 97 | ! determine next virtual temperatures |
|---|
| 98 | tv2 = tconv(kz)*(1.+0.608*qconv(kz)) |
|---|
| 99 | ! interpolate virtual temperature to half-level |
|---|
| 100 | tv = tv1 + (tv2-tv1)*(pconv(kz-1)-phconv(kz))/ & |
|---|
| 101 | (pconv(kz-1)-pconv(kz)) |
|---|
| 102 | if (abs(tv-tvold).gt.0.2) then |
|---|
| 103 | uvzlev(kz) = uvzlev(kz-1) + & |
|---|
| 104 | const*log(pold/pint)* & |
|---|
| 105 | (tv-tvold)/log(tv/tvold) |
|---|
| 106 | else |
|---|
| 107 | uvzlev(kz) = uvzlev(kz-1)+ & |
|---|
| 108 | const*log(pold/pint)*tv |
|---|
| 109 | endif |
|---|
| 110 | tvold=tv |
|---|
| 111 | tv1=tv2 |
|---|
| 112 | pold=pint |
|---|
| 113 | |
|---|
| 114 | end do |
|---|
| 115 | |
|---|
| 116 | |
|---|
| 117 | ktop = 2 |
|---|
| 118 | |
|---|
| 119 | endif ! (if ktop .le. 1) then |
|---|
| 120 | |
|---|
| 121 | ! determine vertical grid position of particle in the eta system |
|---|
| 122 | !**************************************************************** |
|---|
| 123 | |
|---|
| 124 | ztold = ztra1(abs(ipart)) |
|---|
| 125 | ! find old particle grid position |
|---|
| 126 | do kz = 2, nconvtop |
|---|
| 127 | if (uvzlev(kz) .ge. ztold ) then |
|---|
| 128 | levold = kz-1 |
|---|
| 129 | goto 30 |
|---|
| 130 | endif |
|---|
| 131 | end do |
|---|
| 132 | |
|---|
| 133 | ! Particle is above the potentially convective domain. Skip it. |
|---|
| 134 | goto 90 |
|---|
| 135 | |
|---|
| 136 | 30 continue |
|---|
| 137 | |
|---|
| 138 | ! now redistribute particles |
|---|
| 139 | !**************************** |
|---|
| 140 | |
|---|
| 141 | ! Choose a random number and find corresponding level of destination |
|---|
| 142 | ! Random numbers to be evenly distributed in [0,1] |
|---|
| 143 | |
|---|
| 144 | rn = ran3_conv(ipart) |
|---|
| 145 | |
|---|
| 146 | ! initialize levnew |
|---|
| 147 | |
|---|
| 148 | levnew = levold |
|---|
| 149 | |
|---|
| 150 | ffraction = 0. |
|---|
| 151 | totlevmass=dpr(levold)/ga |
|---|
| 152 | do k = 1,nconvtop |
|---|
| 153 | ! for backward runs use the transposed matrix |
|---|
| 154 | if (ldirect.eq.1) then |
|---|
| 155 | ffraction=ffraction+fmassfrac(levold,k) & |
|---|
| 156 | /totlevmass |
|---|
| 157 | else |
|---|
| 158 | ffraction=ffraction+fmassfrac(k,levold) & |
|---|
| 159 | /totlevmass |
|---|
| 160 | endif |
|---|
| 161 | if (rn.le.ffraction) then |
|---|
| 162 | levnew=k |
|---|
| 163 | ! avoid division by zero or a too small number |
|---|
| 164 | ! if division by zero or a too small number happens the |
|---|
| 165 | ! particle is assigned to the center of the grid cell |
|---|
| 166 | if (ffraction.gt.1.e-20) then |
|---|
| 167 | if (ldirect.eq.1) then |
|---|
| 168 | dlevfrac = (ffraction-rn) / fmassfrac(levold,k) * totlevmass |
|---|
| 169 | else |
|---|
| 170 | dlevfrac = (ffraction-rn) / fmassfrac(k,levold) * totlevmass |
|---|
| 171 | endif |
|---|
| 172 | else |
|---|
| 173 | dlevfrac = 0.5 |
|---|
| 174 | endif |
|---|
| 175 | goto 40 |
|---|
| 176 | endif |
|---|
| 177 | end do |
|---|
| 178 | |
|---|
| 179 | 40 continue |
|---|
| 180 | |
|---|
| 181 | ! now assign new position to particle |
|---|
| 182 | |
|---|
| 183 | if (levnew.le.nconvtop) then |
|---|
| 184 | if (levnew.eq.levold) then |
|---|
| 185 | ztra1(abs(ipart)) = ztold |
|---|
| 186 | else |
|---|
| 187 | dlogp = (1.-dlevfrac)* & |
|---|
| 188 | (log(phconv(levnew+1))-log(phconv(levnew))) |
|---|
| 189 | pint = log(phconv(levnew))+dlogp |
|---|
| 190 | dz1 = pint - log(phconv(levnew)) |
|---|
| 191 | dz2 = log(phconv(levnew+1)) - pint |
|---|
| 192 | dz = dz1 + dz2 |
|---|
| 193 | ztra1(abs(ipart)) = (uvzlev(levnew)*dz2+uvzlev(levnew+1)*dz1)/dz |
|---|
| 194 | if (ztra1(abs(ipart)).lt.0.) & |
|---|
| 195 | ztra1(abs(ipart))=-1.*ztra1(abs(ipart)) |
|---|
| 196 | if (ipconv.gt.0) ipconv=-1 |
|---|
| 197 | endif |
|---|
| 198 | endif |
|---|
| 199 | |
|---|
| 200 | ! displace particle according to compensating subsidence |
|---|
| 201 | ! this is done to those particles, that were not redistributed |
|---|
| 202 | ! by the matrix |
|---|
| 203 | !************************************************************** |
|---|
| 204 | |
|---|
| 205 | if (levnew.le.nconvtop.and.levnew.eq.levold) then |
|---|
| 206 | |
|---|
| 207 | ztold = ztra1(abs(ipart)) |
|---|
| 208 | |
|---|
| 209 | ! determine compensating vertical velocity at the levels |
|---|
| 210 | ! above and below the particel position |
|---|
| 211 | ! increase compensating subsidence by the fraction that |
|---|
| 212 | ! is displaced by convection to this level |
|---|
| 213 | |
|---|
| 214 | if (levold.gt.1) then |
|---|
| 215 | temp_levold = tconv(levold-1) + & |
|---|
| 216 | (tconv(levold)-tconv(levold-1)) & |
|---|
| 217 | *(pconv(levold-1)-phconv(levold))/ & |
|---|
| 218 | (pconv(levold-1)-pconv(levold)) |
|---|
| 219 | sub_levold = sub(levold)/(1.-sub(levold)/dpr(levold)*ga) |
|---|
| 220 | wsub(levold)=-1.*sub_levold*r_air*temp_levold/(phconv(levold)) |
|---|
| 221 | else |
|---|
| 222 | wsub(levold)=0. |
|---|
| 223 | endif |
|---|
| 224 | |
|---|
| 225 | temp_levold1 = tconv(levold) + & |
|---|
| 226 | (tconv(levold+1)-tconv(levold)) & |
|---|
| 227 | *(pconv(levold)-phconv(levold+1))/ & |
|---|
| 228 | (pconv(levold)-pconv(levold+1)) |
|---|
| 229 | sub_levold1 = sub(levold+1)/(1.-sub(levold+1)/dpr(levold+1)*ga) |
|---|
| 230 | wsub(levold+1)=-1.*sub_levold1*r_air*temp_levold1/ & |
|---|
| 231 | (phconv(levold+1)) |
|---|
| 232 | |
|---|
| 233 | ! interpolate wsub to the vertical particle position |
|---|
| 234 | |
|---|
| 235 | dz1 = ztold - uvzlev(levold) |
|---|
| 236 | dz2 = uvzlev(levold+1) - ztold |
|---|
| 237 | dz = dz1 + dz2 |
|---|
| 238 | |
|---|
| 239 | wsubpart = (dz2*wsub(levold)+dz1*wsub(levold+1))/dz |
|---|
| 240 | ztra1(abs(ipart)) = ztold+wsubpart*real(lsynctime) |
|---|
| 241 | if (ztra1(abs(ipart)).lt.0.) then |
|---|
| 242 | ztra1(abs(ipart))=-1.*ztra1(abs(ipart)) |
|---|
| 243 | endif |
|---|
| 244 | |
|---|
| 245 | endif !(levnew.le.nconvtop.and.levnew.eq.levold) |
|---|
| 246 | |
|---|
| 247 | ! Maximum altitude .5 meter below uppermost model level |
|---|
| 248 | !******************************************************* |
|---|
| 249 | |
|---|
| 250 | 90 continue |
|---|
| 251 | |
|---|
| 252 | if (ztra1(abs(ipart)) .gt. height(nz)-0.5) & |
|---|
| 253 | ztra1(abs(ipart)) = height(nz)-0.5 |
|---|
| 254 | |
|---|
| 255 | end subroutine redist |
|---|
