Changeset 3c1da52 in flexpart.git for src/wetdepo.f90
- Timestamp:
- Jul 15, 2015, 9:45:31 AM (9 years ago)
- Branches:
- svn-petra
- Children:
- c07b09e
- Parents:
- 31674b5
- git-author:
- Ignacio Pisso <Ignacio.Pisso@…> (07/15/15 09:23:56)
- git-committer:
- Ignacio Pisso <Ignacio.Pisso@…> (07/15/15 09:45:31)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
src/wetdepo.f90
r4fbe7a5 r3c1da52 1 1 !********************************************************************** 2 ! Copyright 1998 ,1999,2000,2001,2002,2005,2007,2008,2009,2010*2 ! Copyright 1998-2015 * 3 3 ! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, * 4 4 ! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann * … … 39 39 ! Code may not be correct for decay of deposition! * 40 40 ! * 41 ! PS, 2/2015: implement wet depo quick fix from 2011/12 in this version 42 ! it implements interpolated and improved clouds 43 ! Also, certain deficiencies for thin clouds fixed also 44 ! Pass itage to wetdepokernel 45 ! * 41 46 !***************************************************************************** 42 47 ! * … … 71 76 72 77 integer :: jpart,itime,ltsample,loutnext,ldeltat,i,j,ix,jy 73 integer :: ngrid,itage,nage, hz,il,interp_time, n, clouds_v74 integer :: ks, kp 78 integer :: ngrid,itage,nage,kz,il,interp_time,n 79 integer :: ks, kp, n1,n2, icbot,ictop, indcloud 75 80 real :: S_i, act_temp, cl, cle ! in cloud scavenging 76 real :: clouds_h ! cloud height for the specific grid point 77 real :: xtn,ytn,lsp,convp,cc,grfraction,prec,wetscav 81 real :: xtn,ytn,lsp,convp,cc,grfraction,prec,wetscav,wetscavold,precsub,f 78 82 real :: wetdeposit(maxspec),restmass 79 83 real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled … … 97 101 !************************ 98 102 103 particle_loop: & 99 104 do jpart=1,numpart 100 105 101 106 if (itra1(jpart).eq.-999999999) goto 20 102 if (ldirect.eq.1)then107 if (ldirect.eq.1) then 103 108 if (itra1(jpart).gt.itime) goto 20 104 109 else 105 110 if (itra1(jpart).lt.itime) goto 20 106 111 endif 112 107 113 ! Determine age class of the particle 108 114 itage=abs(itra1(jpart)-itramem(jpart)) … … 123 129 goto 23 124 130 endif 125 end 131 enddo 126 132 23 continue 127 133 … … 131 137 132 138 if (ngrid.gt.0) then 133 xtn= (xtra1(jpart)-xln(ngrid))*xresoln(ngrid)134 ytn= (ytra1(jpart)-yln(ngrid))*yresoln(ngrid)139 xtn=real(xtra1(jpart)-xln(ngrid))*xresoln(ngrid) 140 ytn=real(ytra1(jpart)-yln(ngrid))*yresoln(ngrid) 135 141 ix=int(xtn) 136 142 jy=int(ytn) … … 148 154 149 155 if (ngrid.eq.0) then 150 call interpol_rain(lsprec,convprec,tcc,nxmax,nymax, 151 1,nx,ny,memind,real(xtra1(jpart)),real(ytra1(jpart)),1, &152 memtime(1),memtime(2),interp_time,lsp,convp,cc)156 call interpol_rain(lsprec,convprec,tcc,nxmax,nymax,1,nx,ny,memind,& 157 real(xtra1(jpart)),real(ytra1(jpart)), & 158 1,memtime(1),memtime(2),interp_time,lsp,convp,cc,icbot,ictop) 153 159 else 154 160 call interpol_rain_nests(lsprecn,convprecn,tccn, & 155 nxmaxn,nymaxn,1,maxnests,ngrid,nxn,nyn,memind,xtn,ytn,1, & 156 memtime(1),memtime(2),interp_time,lsp,convp,cc) 157 endif 158 159 if ((lsp.lt.0.01).and.(convp.lt.0.01)) goto 20 161 nxmaxn,nymaxn,1,maxnests,ngrid,nxn,nyn,memind,xtn,ytn, & 162 1,memtime(1),memtime(2),interp_time,lsp,convp,cc,icbot,ictop) 163 endif 164 165 ! PS 2012/2015: subtract a small value, eg 0.01 mm/h, 166 ! to remove spurious precip; replaces previous code 167 prec = lsp+convp 168 precsub = 0.01 169 if (prec .lt. precsub) then 170 goto 20 171 else 172 f = (prec-precsub)/prec 173 lsp = f*lsp 174 convp = f*convp 175 endif 176 160 177 ! get the level were the actual particle is in 161 178 do il=2,nz 162 179 if (height(il).gt.ztra1(jpart)) then 163 hz=il-1180 kz=il-1 164 181 goto 26 165 182 endif 166 end 183 enddo 167 184 26 continue 168 185 169 186 n=memind(2) 170 187 if (abs(memtime(1)-interp_time).lt.abs(memtime(2)-interp_time)) & 171 n=memind(1)188 n=memind(1) 172 189 173 190 ! if there is no precipitation or the particle is above the clouds no 174 191 ! scavenging is done 175 176 if (ngrid.eq.0) then 177 clouds_v=clouds(ix,jy,hz,n) 178 clouds_h=cloudsh(ix,jy,n) 179 else 180 clouds_v=cloudsn(ix,jy,hz,n,ngrid) 181 clouds_h=cloudsnh(ix,jy,n,ngrid) 182 endif 183 !write(*,*) 'there is 184 ! + precipitation',(clouds(ix,jy,ihz,n),ihz=1,20),lsp,convp,hz 185 if (clouds_v.le.1) goto 20 186 !write (*,*) 'there is scavenging' 192 ! PS: part of 2011/2012/2015 fix, replaces previous code 193 194 if (ztra1(jpart) .le. float(ictop)) then 195 if (ztra1(jpart) .gt. float(icbot)) then 196 indcloud = 2 ! in-cloud 197 else 198 indcloud = 1 ! below-cloud 199 endif 200 elseif (ictop .eq. icmv) then 201 indcloud = 0 ! no cloud found, use old scheme 202 else 203 goto 20 ! above cloud 204 endif 205 187 206 188 207 ! 1) Parameterization of the the area fraction of the grid cell where the … … 228 247 !********************************************************** 229 248 249 species_loop: & 230 250 do ks=1,nspec ! loop over species 251 231 252 wetdeposit(ks)=0. 232 253 wetscav=0. 233 254 234 255 ! NIK 09.12.13: allowed to turn off either below cloud or in-cloud by including TWO separate blocks of if tests 235 236 237 !****i******************* 238 ! BELOW CLOUD SCAVENGING 239 !*********************** 240 if (weta(ks).gt.0. .and. clouds_v.ge.4) then !if positive below-cloud coefficient given in SPECIES file, and if in below cloud regime 241 ! for aerosols and not highliy soluble substances weta=5E-6 242 wetscav=weta(ks)*prec**wetb(ks) ! scavenging coefficient 243 endif !end below-cloud 244 245 246 !******************** 247 ! IN CLOUD SCAVENGING 248 !******************** 249 if (weta_in(ks).gt.0. .and. clouds_v.lt.4) then !if positive in-cloud coefficient given in SPECIES file, and if in in-cloud regime 250 251 if (ngrid.gt.0) then 252 act_temp=ttn(ix,jy,hz,n,ngrid) 253 else 254 act_temp=tt(ix,jy,hz,n) 255 endif 256 if (weta(ks).gt.0.) then 257 ! positive below-cloud coefficient from SPECIES file 258 259 if (indcloud .eq. 1) then ! below-cloud scavenging 260 261 ! for aerosols and not highly soluble substances weta=5E-6 262 wetscav=weta(ks)*prec**wetb(ks) ! scavenging coefficient 263 264 elseif (indcloud .eq. 2) then ! in-cloud scavenging 265 266 if (ngrid.gt.0) then 267 act_temp=ttn(ix,jy,kz,n,ngrid) 268 else 269 act_temp=tt(ix,jy,kz,n) 270 endif 256 271 257 272 ! NIK 31.01.2013: SPECIES defined parameters for the in-cloud scavening … … 259 274 ! wetb_in=0.36 (default) 260 275 ! wetc_in=0.9 (default) 261 ! wetd_in: Scaling factor for the total in-cloud scavenging (default 1.0 -no scaling)276 ! wetd_in: Scaling factor for the total in-cloud scavenging (default 1.0 - no scaling) 262 277 cl=weta_in(ks)*prec**wetb_in(ks) 263 278 if (dquer(ks).gt.0) then ! is particle 264 279 S_i=wetc_in(ks)/cl 265 280 else ! is gas 266 cle=(1 -cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl267 S_i=1 /cle281 cle=(1.-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl 282 S_i=1./cle 268 283 endif 269 wetscav=S_i*prec/3.6E6/clouds_h/wetd_in(ks) 270 ! write(*,*) 'in. wetscav:' 271 ! + ,wetscav,cle,cl,act_temp,prec,clouds_h 272 273 endif ! end in-cloud 274 275 284 wetscav=S_i*prec/3.6E6/wetd_in(ks)/(ictop-icbot) ! 3.6e6 converts mm/h to m/s 285 ! PS - prevent extrem values of wetscav: 286 wetscavold = 2.e-5*prec**0.8 287 wetscav = min( 0.1*wetscav, wetscavold ) ! here we sacle the above wetscav. This is now twice, also in wetd_in. 288 289 else ! PS: no cloud diagnosed, old scheme, 290 291 ! PS using with fixed a,b for simplicity, one may wish to change!! 292 wetscav = 2.e-5*prec**0.8 ! was before: 1.e-4*prec**0.62 293 294 endif ! end regime cases 276 295 277 !************************************************** 278 ! CALCULATE DEPOSITION 279 !************************************************** 280 ! if (wetscav.le.0) write (*,*) 'neg, or 0 wetscav!' 281 ! + ,wetscav,cle,cl,act_temp,prec,clouds_h,clouds_v 282 283 if (wetscav.gt.0.) then 284 wetdeposit(ks)=xmass1(jpart,ks)* & 285 (1.-exp(-wetscav*abs(ltsample)))*grfraction ! wet deposition 286 else ! if no scavenging 287 wetdeposit(ks)=0. 288 endif 289 290 296 ! calculate deposition 297 wetdeposit(ks)=xmass1(jpart,ks)* & 298 (1.-exp(-wetscav*abs(ltsample)))*grfraction 291 299 restmass = xmass1(jpart,ks)-wetdeposit(ks) 292 300 if (ioutputforeachrelease.eq.1) then … … 297 305 if (restmass .gt. smallnum) then 298 306 xmass1(jpart,ks)=restmass 299 ! depostatistic 300 ! wetdepo_sum(ks,kp)=wetdepo_sum(ks,kp)+wetdeposit(ks) 301 ! depostatistic 307 ! depostatistic: 308 !! wetdepo_sum(ks,kp)=wetdepo_sum(ks,kp)+wetdeposit(ks) 302 309 else 303 310 xmass1(jpart,ks)=0. … … 305 312 ! Correct deposited mass to the last time step when radioactive decay of 306 313 ! gridded deposited mass was calculated 307 if (decay(ks).gt.0.) then314 if (decay(ks).gt.0.) & 308 315 wetdeposit(ks)=wetdeposit(ks)*exp(abs(ldeltat)*decay(ks)) 309 endif 310 311 312 313 end do !all species 314 315 ! Sabine Eckhardt, June 2008 create deposition runs only for forward runs 316 ! Add the wet deposition to accumulated amount on output grid and nested output grid 317 !***************************************************************************** 316 317 else ! weta(k) <= 0 318 319 wetdeposit(ks)=0. 320 321 endif 322 323 enddo species_loop 324 325 ! Sabine Eckhardt, June 2008: write deposition only in forward runs 326 ! add the wet deposition from this step to accumulated amount 327 ! on output grid and nested output grid 318 328 319 329 if (ldirect.eq.1) then 320 call wetdepokernel(nclass(jpart),wetdeposit, real(xtra1(jpart)),&321 real(ytra1(jpart)),nage,kp)330 call wetdepokernel(nclass(jpart),wetdeposit, & 331 real(xtra1(jpart)),real(ytra1(jpart)),itage,nage,kp) 322 332 if (nested_output.eq.1) call wetdepokernel_nest(nclass(jpart), & 323 wetdeposit,real(xtra1(jpart)),real(ytra1(jpart)), nage,kp)324 endif 325 326 20 continue 327 end do ! all particles333 wetdeposit,real(xtra1(jpart)),real(ytra1(jpart)),itage,nage,kp) 334 endif 335 336 20 continue ! jump here for particles not to be treated 337 enddo particle_loop 328 338 329 339 end subroutine wetdepo
Note: See TracChangeset
for help on using the changeset viewer.