source: branches/jerome/src_flexwrf_v3.1/advance.f90 @ 16

Last change on this file since 16 was 16, checked in by jebri, 10 years ago

sources for flexwrf v3.1

File size: 51.9 KB
Line 
1!***********************************************************************
2!* Copyright 2012,2013                                                *
3!* Jerome Brioude, Delia Arnold, Andreas Stohl, Wayne Angevine,       *
4!* John Burkhart, Massimo Cassiani, Adam Dingwell, Richard C Easter, Sabine Eckhardt,*
5!* Stephanie Evan, Jerome D Fast, Don Morton, Ignacio Pisso,          *
6!* Petra Seibert, Gerard Wotawa, Caroline Forster, Harald Sodemann,   *
7!* This file is part of FLEXPART WRF                                   *
8!*                                                                     *
9!* FLEXPART is free software: you can redistribute it and/or modify    *
10!* it under the terms of the GNU General Public License as published by*
11!* the Free Software Foundation, either version 3 of the License, or   *
12!* (at your option) any later version.                                 *
13!*                                                                     *
14!* FLEXPART is distributed in the hope that it will be useful,         *
15!* but WITHOUT ANY WARRANTY; without even the implied warranty of      *
16!* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
17!* GNU General Public License for more details.                        *
18!*                                                                     *
19!* You should have received a copy of the GNU General Public License   *
20!* along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
21!***********************************************************************
22
23      subroutine advance(itime,nrelpoint,ldt,up,vp,wp, &
24      usigold,vsigold,wsigold,nstop,xt,yt,zt,prob,icbt, &
25      ngrid,depoindicator,indzindicator,cpt2,ompid,myid,n_threads,mts)  !comment by mc: ...,ompid,myid,n_threads) added  n_threads for MT parallel random number generator
26
27!                          i    i  i/oi/oi/o
28!       i/o     i/o     i/o     o  i/oi/oi/o i/o  i/o
29!*******************************************************************************
30!                                                                              *
31!  Note:  This is the FLEXPART_WRF version of subroutine gridcheck.            *
32!    The computational grid is the WRF x-y grid rather than lat-lon.           *
33!                                                                              *
34!  Calculation of turbulent particle trajectories utilizing a                  *
35!  zero-acceleration scheme, which is corrected by a numerically more accurate *
36!  Petterssen scheme whenever possible.                                        *
37!                                                                              *
38!  Particle positions are read in, incremented, and returned to the calling    *
39!  program.                                                                    *
40!                                                                              *
41!  In different regions of the atmosphere (PBL vs. free troposphere),          *
42!  different parameters are needed for advection, parameterizing turbulent     *
43!  velocities, etc. For efficiency, different interpolation routines have      *
44!  been written for these different cases, with the disadvantage that there    *
45!  exist several routines doing almost the same. They all share the            *
46!  included file 'includeinterpol'. The following                              *
47!  interpolation routines are used:                                            *
48!                                                                              *
49!  interpol_all(_nests)     interpolates everything (called inside the PBL)    *
50!  interpol_misslev(_nests) if a particle moves vertically in the PBL,         *
51!                           additional parameters are interpolated if it       *
52!                           crosses a model level                              *
53!  interpol_wind(_nests)    interpolates the wind and determines the           *
54!                           standard deviation of the wind (called outside PBL)*
55!                           also interpolates potential vorticity              *
56!  interpol_wind_short(_nests) only interpolates the wind (needed for the      *
57!                           Petterssen scheme)                                 *
58!  interpol_vdep(_nests)    interpolates deposition velocities                 *
59!                                                                              *
60!                                                                              *
61!     Author: A. Stohl                                                         *
62!                                                                              *
63!     16 December 1997                                                         *
64!                                                                              *
65!  Changes:                                                                    *
66!                                                                              *
67!  8 April 2000: Deep convection parameterization                              *
68!                                                                              *
69!  May 2002: Petterssen scheme introduced                                      *
70!                                                                              *
71!  26 Oct 2005, R. Easter - changes for horizontal grid in m instead of lat,lon*
72!  10 Nov 2005, R. Easter - zero turbulent wind components is                  *
73!                           turbulence is turned off                           *
74!  Mar 2012, J. Brioude: modification to handle openmp.                        *
75!  turbulence option 3 is not going   *
76!  to work. but it shouldn't be used anyway ^M
77!  Jan 2013  M. Cassiani (look for mc or MC in the code):
78!  *^M
79!               introduction of CBL skewed turbulence model
80!               *^M
81!               & parallel random number generation
82!               *                       *
83!********************************************************************************
84!*******************************************************************************
85!                                                                              *
86! Variables:                                                                   *
87! cbt                1 if particle not transferred to forbidden state, else -1 *
88! dawsave            accumulated displacement in along-wind direction          *
89! dcwsave            accumulated displacement in cross-wind direction          *
90! dxsave             accumulated displacement in longitude                     *
91! dysave             accumulated displacement in latitude                      *
92! h [m]              Mixing height                                             *
93! lwindinterv [s]    time interval between two wind fields                     *
94! itime [s]          time at which this subroutine is entered                  *
95! itimec [s]         actual time, which is incremented in this subroutine      *
96! href [m]           height for which dry deposition velocity is calculated    *
97! ladvance [s]       Total integration time period                             *
98! ldirect            1 forward, -1 backward                                    *
99! ldt [s]            Time step for the next integration                        *
100! lsynctime [s]      Synchronisation interval of FLEXPART                      *
101! ngrid              index which grid is to be used                            *
102! nrand              index for a variable to be picked from rannumb            *
103! nstop              if > 1 particle has left domain and must be stopped       *
104! prob               probability of absorption due to dry deposition           *
105! rannumb(maxrand)   normally distributed random variables                     *
106! rhoa               air density                                               *
107! rhograd            vertical gradient of the air density                      *
108! up,vp,wp           random velocities due to turbulence (along wind, cross    *
109!                    wind, vertical wind                                       *
110! usig,vsig,wsig     mesoscale wind fluctuations                               *
111! usigold,vsigold,wsigold  like usig, etc., but for the last time step         *
112! vdepo              Deposition velocities for all species                     *
113! xt,yt,zt           Particle position                                         *
114!                                                                              *
115!*******************************************************************************
116
117  use point_mod
118  use par_mod
119  use com_mod
120  use mt_stream  !added by mc for random number generation^M
121!  use test_well_mod !added by mc for testting well mixed
122!  use interpol_mod
123!  use hanna_mod
124  use cmapf_mod
125!  use ieee_arithmetic
126! include 'sprng_f.h'
127!  use ran_mod
128!      include 'includepar'
129!      include 'includecom'
130!      include 'includeinterpol'
131!      include 'includehanna'
132
133  implicit none
134  real(kind=dp) :: xt,yt
135  real :: zt,xts,yts,weight
136  integer :: itime,itimec,nstop,ldt,i,j,k,nrand,loop,memindnext
137  integer :: ngr,nix,njy,ks,nsp,nrelpoint,ii,ompid,myid,nombre
138  real :: dz,dz1,dz2,xlon,ylat,xpol,ypol,gridsize
139  real :: ru,rv,rw,dt,ux,vy,cosfact,xtn,ytn,tropop
140  real :: prob(maxspec),up,vp,wp,dxsave,dysave,dawsave
141  real :: dcwsave,mu,mv
142  real :: usigold,vsigold,wsigold,r,rs
143  real :: uold,vold,wold,vdepo(maxspec)
144  !real uprof(nzmax),vprof(nzmax),wprof(nzmax)
145  !real usigprof(nzmax),vsigprof(nzmax),wsigprof(nzmax)
146  !real rhoprof(nzmax),rhogradprof(nzmax)
147  real :: rhoa,rhograd,ran3,delz,dtf,rhoaux,dtftlw,uxscale,wpscale
148  real :: ptot_lhh,Q_lhh,phi_lhh,ath,bth   !modified by mc added for CBL scheme
149  real :: old_wp_buf,del_test              !modified by mc added for CBL scheme re-initlization fo particle after NaN
150  integer(kind=2) :: icbt
151  real,parameter :: eps=nxmax/3.e5,eps2=1.e-9
152  integer :: flagrein                      !re-initialization flag for particles: modified by mc
153
154  real :: uprof(nzmax),vprof(nzmax),wprof(nzmax)
155  real :: usigprof(nzmax),vsigprof(nzmax),wsigprof(nzmax)
156  real :: rhoprof(nzmax),rhogradprof(nzmax)
157  real :: tkeprof(nzmax),pttprof(nzmax)
158  real :: u,v,w,usig,vsig,wsig,pvi
159  real(kind=dp) :: xtold
160  real :: p1,p2,p3,p4,ddx,ddy,rddx,rddy,dtt,dt1,dt2
161  integer :: ix,jy,ixp,jyp,ngrid,indz,indzp,cpt2,maxrand2
162  logical :: depoindicator(maxspec)
163  logical :: indzindicator(nzmax)
164
165  real :: ust,wst,ol,h,zeta,sigu,sigv,tlu,tlv,tlw
166  real :: sigw,dsigwdz,dsigw2dz
167
168  real :: wp2, zt2, ust2, wst2, h2, rhoa2, rhograd2, sigw2, &
169    dsigwdz2, tlw2, ptot_lhh2, Q_lhh2, phi_lhh2, ath2, bth2, ol2
170
171  logical :: isnan2
172
173!!! CHANGE: TEST OF THE WELL-MIXED CRITERION
174!        integer iclass
175!        parameter(iclass=10)
176!        double precision zacc,tacc,t(iclass),th(0:iclass),hsave
177!        logical dump
178!        save zacc,tacc,t,th,hsave,dump
179!c        itimeod=0
180!!! CHANGE
181
182  integer :: idummy = 7
183  real    :: settling = 0.
184  !added by mc for random number generation --------------------- 
185  integer ::  n_threads !added by mc for parallel random number generation
186  integer(4) :: rannum
187  real(4) :: real_rannum
188  type (mt_state) :: mts (0: MAX_STREAM)
189  integer,SAVE :: nan_count(max_STREAM)=0
190!-------------------------------------------------------
191
192
193!!! CHANGE: TEST OF THE WELL-MIXED CRITERION
194!      if (idummy.eq.-7) then
195!      open(550,file='WELLMIXEDTEST')
196!      do 17 i=0,iclass
197!17      th(i)=real(i)/real(iclass)
198!      endif
199!!! CHANGE
200
201!    if (nombre.eq.103) print*,'usig -1',usig,xts,yts,zt
202     if (xt.ne.xt .or. abs(xt).gt.1000.) print*,'problem 0', xt,yt,zt,itime,myid,ompid
203     xtold=xt
204!   print *,'aa',xt,yt,zt,u,v,w
205      nstop=0
206      do i=1,nmixz
207        indzindicator(i)=.true.
208       enddo
209
210      if (DRYDEP) then    ! reset probability for deposition
211        do  ks=1,nspec
212          depoindicator(ks)=.true.
213          prob(ks)=0.
214       enddo
215      endif
216
217      dxsave=0.           ! reset position displacements
218      dysave=0.           ! due to mean wind
219      dawsave=0.          ! and turbulent wind
220      dcwsave=0.
221      usig=0.
222      vsig=0.
223      wsig=0.
224       ust=0.
225      wst=0.
226      ol=0.
227      h=0.
228      zeta=0.
229      sigu=0.
230      sigv=0.
231      tlu=0.
232      tlv=0.
233      tlw=0.
234      sigw=0.
235      dsigwdz=0.
236      dsigw2dz=0.
237!     wp=0.
238        itimec=itime
239      idummy=7
240       if (newrandomgen.eq.0) then
241!      cpt2=cpt2+ompid+1
242      cpt2=cpt2+1
243!     cpt2=cpt2+1000367
244      cpt2=mod(cpt2,maxrand)+1;
245
246!      nrand=int(ran3(idummy,inext,inextp,ma,iff)*real(maxrand-1))+1
247!     nrand=cpt
248      nrand=cpt2+ompid*maxrand
249!     print*,cpt2,maxrand,maxrandomp,maxomp
250!       print*, rannumb(nrand),nrelpoint
251!     print*,rannumb(nrand),myid,ompid
252!      if (nrelpoint.ge.993 .and. nrelpoint.le.998) then
253!      write(22,*),itime,cpt2,rannumb(nrand),nrelpoint
254!!,myid,OMP_GET_THREAD_NUM()
255 !!      write(22,*),itime,cpt2,rannumb(nrand),nrelpoint
256!         endif
257
258
259!        if (nrand+2.gt.maxrand) nrand=1
260!    print*,rannumb(nrand)
261       maxrand2=maxrandomp
262
263        else
264!-------------------------------------------------------------------------------------------------
265!----- added by MC: parallel random nuymber generation using MT generator ------------------------
266       !print *,'varie3',ompid,myid,n_threads
267!      rannum=genrand_int32(mts(ompid+1+(myid*n_threads)))  !integer random number at 32 bit resolution
268       rannum=genrand_int32(mts(ompid+1))  !integer random number at 32 bit resolution
269       real_rannum = sngl(0.5_DP + 0.2328306e-9_DP * rannum) !conversion to single precision 32bit real
270       nrand=int(real_rannum*real(maxrand-1))+1
271        !print *,'varie4',rannum,real_rannum,nrand
272!--------------------------------------------------------------------------------------------------
273        maxrand2=maxrand
274
275        end if
276       
277! Determine whether lat/long grid or polarstereographic projection
278! is to be used
279! Furthermore, determine which nesting level to be used
280!*****************************************************************
281
282      if (nglobal.and.(yt.gt.switchnorthg)) then
283        write(*,*)
284        write(*,*) '*** stopping in advance ***'
285        write(*,*) '    the n-pole code section should not be active'
286        write(*,*)
287        ngrid=-1
288      else if (sglobal.and.(yt.lt.switchsouthg)) then
289        write(*,*)
290        write(*,*) '*** stopping in advance ***'
291        write(*,*) '    the s-pole code section should not be active'
292        write(*,*)
293        ngrid=-2
294      else
295        ngrid=0
296        do j=numbnests,1,-1
297          if ((xt.gt.xln(j)+eps).and.(xt.lt.xrn(j)-eps).and. &
298          (yt.gt.yln(j)+eps).and.(yt.lt.yrn(j)-eps)) then
299            ngrid=j
300            goto 23
301          endif
302        enddo
30323      continue
304      endif
305
306
307!***************************
308! Interpolate necessary data
309!***************************
310
311      if (abs(itime-memtime(1)).lt.abs(itime-memtime(2))) then
312        memindnext=1
313      else
314        memindnext=2
315      endif
316
317! Determine nested grid coordinates
318!**********************************
319
320      if (ngrid.gt.0) then
321        xtn=(xt-xln(ngrid))*xresoln(ngrid)
322        ytn=(yt-yln(ngrid))*yresoln(ngrid)
323        ix=int(xtn)
324        jy=int(ytn)
325    nix=nint(xtn)
326    njy=nint(ytn)
327
328      else
329        ix=int(xt)
330        jy=int(yt)
331    nix=nint(xt)
332    njy=nint(yt)
333      endif
334      ixp=ix+1
335      jyp=jy+1
336
337     if (ix.lt.0) print*,'problem', xt,xtold,yt,zt,itime,myid,ompid,nrelpoint
338
339! Compute maximum mixing height around particle position
340!*******************************************************
341
342      h=0.
343      if (ngrid.le.0) then
344        do k=1,2
345          do j=jy,jyp
346            do i=ix,ixp
347              if (hmix(i,j,1,k).gt.h) h=hmix(i,j,1,k)
348        end do
349      end do
350    end do
351    tropop=tropopause(nix,njy,1,1)
352  else
353    do k=1,2
354      do j=jy,jyp
355        do i=ix,ixp
356          if (hmixn(i,j,1,k,ngrid).gt.h) h=hmixn(i,j,1,k,ngrid)
357        end do
358      end do
359    end do
360    tropop=tropopausen(nix,njy,1,1,ngrid)
361  endif
362
363  zeta=zt/h
364
365
366!*************************************************************
367! If particle is in the PBL, interpolate once and then make a
368! time loop until end of interval is reached
369!*************************************************************
370!      print*,'zeta',zeta,h,zt,xt
371      if (zeta.le.1.) then
372
373! BEGIN TIME LOOP
374!================
375
376        loop=0
377100       loop=loop+1
378          if (method.eq.1) then
379            ldt=min(ldt,abs(lsynctime-itimec+itime))
380            itimec=itimec+ldt*ldirect
381          else
382            ldt=abs(lsynctime)
383            itimec=itime+lsynctime
384          endif
385          dt=real(ldt)
386
387          zeta=zt/h
388
389
390!   print *,'xx0',OMP_GET_THREAD_NUM(),loop,xt,yt,zt,xts,yts,u,v,w
391          if (loop.eq.1) then
392!    if (nombre.eq.103) print*,'usig 0',usig,xt,yt,zt
393            if (ngrid.le.0) then
394              xts=real(xt)
395              yts=real(yt)
396!    if (nombre.eq.103) print*,'usig 0',usig,xts,yts,zt
397              call interpol_all(itime,xts,yts,zt, &
398    uprof,vprof,wprof, usigprof,vsigprof,wsigprof, &
399    rhoprof,rhogradprof, tkeprof,pttprof, &
400    u,v,w,usig,vsig,wsig,pvi, &
401    p1,p2,p3,p4,ddx,ddy,rddx,rddy,dtt,dt1,dt2, &
402    ix,jy,ixp,jyp,ngrid,indz,indzp, depoindicator, &
403    indzindicator, &
404    ust,wst,ol,h,zeta,sigu,sigv,tlu,tlv,tlw, &
405    sigw,dsigwdz,dsigw2dz,mu,mv)
406
407            else
408              call interpol_all_nests(itime,xtn,ytn,zt, &
409    uprof,vprof,wprof, usigprof,vsigprof,wsigprof, &
410    rhoprof,rhogradprof, tkeprof,pttprof, &
411    u,v,w,usig,vsig,wsig,pvi, &
412    p1,p2,p3,p4,ddx,ddy,rddx,rddy,dtt,dt1,dt2, &
413    ix,jy,ixp,jyp,ngrid,indz,indzp, depoindicator, &
414    indzindicator, &
415    ust,wst,ol,h,zeta,sigu,sigv,tlu,tlv,tlw, &
416    sigw,dsigwdz,dsigw2dz,mu,mv)
417            endif
418!    if (nombre.eq.103) print*,'usig 1',usig,xts,yts,zt
419
420          else
421
422
423!   print *,'xx',OMP_GET_THREAD_NUM(),xt,yt,zt,xts,yts,u,v,w
424! Determine the level below the current position for u,v,rho
425!***********************************************************
426
427            do i=2,nz
428              if (height(i).gt.zt) then
429                indz=i-1
430                indzp=i
431                goto 6
432              endif
433              enddo
4346           continue
435
436! If one of the levels necessary is not yet available,
437! calculate it
438!*****************************************************
439
440            do i=indz,indzp
441              if (indzindicator(i)) then
442                if (ngrid.le.0) then
443!    if (nombre.eq.103) print*,'in usig 2',usig
444                  call interpol_misslev(i,xt,yt,zt, &
445    uprof,vprof,wprof, usigprof,vsigprof,wsigprof, &
446    rhoprof,rhogradprof, tkeprof,pttprof, &
447    u,v,w,usig,vsig,wsig,pvi, &
448    p1,p2,p3,p4,ddx,ddy,rddx,rddy,dtt,dt1,dt2, &
449    ix,jy,ixp,jyp,ngrid,indz,indzp, depoindicator, &
450    indzindicator, &
451    ust,wst,ol,h,zeta,sigu,sigv,tlu,tlv,tlw, &
452    sigw,dsigwdz,dsigw2dz)
453!JB mw not needed here
454                else
455                  call interpol_misslev_nests(i,xt,yt,zt, &
456    uprof,vprof,wprof, usigprof,vsigprof,wsigprof, &
457    rhoprof,rhogradprof, tkeprof,pttprof, &
458    u,v,w,usig,vsig,wsig,pvi, &
459    p1,p2,p3,p4,ddx,ddy,rddx,rddy,dtt,dt1,dt2, &
460    ix,jy,ixp,jyp,ngrid,indz,indzp, depoindicator, &
461    indzindicator, &
462    ust,wst,ol,h,zeta,sigu,sigv,tlu,tlv,tlw, &
463    sigw,dsigwdz,dsigw2dz)
464                endif
465              endif
466            enddo
467          endif
468
469!    if (nombre.eq.103) print*,'usig 2',usig
470
471! Vertical interpolation of u,v,w,rho and drhodz
472!***********************************************
473
474! Vertical distance to the level below and above current position
475! both in terms of (u,v) and (w) fields
476!****************************************************************
477
478          dz=1./(height(indzp)-height(indz))
479          dz1=(zt-height(indz))*dz
480          dz2=(height(indzp)-zt)*dz
481
482          u=dz1*uprof(indzp)+dz2*uprof(indz)
483          v=dz1*vprof(indzp)+dz2*vprof(indz)
484          w=dz1*wprof(indzp)+dz2*wprof(indz)
485          rhoa=dz1*rhoprof(indzp)+dz2*rhoprof(indz)
486          rhograd=dz1*rhogradprof(indzp)+dz2*rhogradprof(indz)
487
488
489! Compute the turbulent disturbances
490! Determine the sigmas and the timescales
491!****************************************
492          if (turb_option .eq. turb_option_mytke) then
493! FLEXPART-WRF
494
495!               write(*,*)'itime=',itime,'xt=',xt,'yt=',yt
496              call tke_partition_my(zt, &
497   ust,wst,ol,h,zeta,sigu,sigv,tlu,tlv,tlw, &
498   sigw,dsigwdz,dsigw2dz,uprof,vprof,tkeprof,pttprof,indz,indzp)
499          else if (turb_option .eq. turb_option_tke) then
500              call tke_partition_hanna(zt, &
501   ust,wst,ol,h,zeta,sigu,sigv,tlu,tlv,tlw, &
502   sigw,dsigwdz,dsigw2dz,uprof,vprof,tkeprof,pttprof,indz,indzp)
503          else
504             if (turbswitch) then
505               call hanna(zt,  &
506    ust,wst,ol,h,zeta,sigu,sigv,tlu,tlv,tlw, &
507    sigw,dsigwdz,dsigw2dz)
508                       
509             else
510               call hanna1(zt, &
511    ust,wst,ol,h,zeta,sigu,sigv,tlu,tlv,tlw, &
512    sigw,dsigwdz,dsigw2dz)
513
514            endif
515          endif   
516!      print*,   ust,wst,ol,h,zeta,sigu,sigv,tlu,tlv,tlw, &
517!  sigw,dsigwdz,dsigw2dz,indz,indzp
518!JB
519!          if (h/abs(ol).lt.1.) then
520!c        print*,itime,'h and ol',h,ol,'neutral'
521!           reflect_switch=0
522!       else if (ol.lt.0.) then
523!c        print*,itime,'h and ol',h,ol,'unstable'
524!c          reflect_switch=1
525!           reflect_switch=0
526!        else
527!c        print*,itime,'h and ol',h,ol,'stable'
528!           reflect_switch=0
529!         endif
530!*****************************************
531! Determine the new diffusivity velocities
532!*****************************************
533
534! Horizontal components
535!**********************
536
537!         if (nrand+1.gt.maxrandomp) nrand=1
538          if (nrand+1.gt.maxrand2) nrand=1
539          if (dt/tlu.lt..5) then
540            up=(1.-dt/tlu)*up+rannumb(nrand)*sigu*sqrt(2.*dt/tlu)
541          else
542            ru=exp(-dt/tlu)
543            up=ru*up+rannumb(nrand)*sigu*sqrt(1.-ru**2)
544          endif
545          if (dt/tlv.lt..5) then
546            vp=(1.-dt/tlv)*vp+rannumb(nrand+1)*sigv*sqrt(2.*dt/tlv)
547          else
548            rv=exp(-dt/tlv)
549            vp=rv*vp+rannumb(nrand+1)*sigv*sqrt(1.-rv**2)
550          endif
551          nrand=nrand+2
552
553
554!         if (nrand+ifine.gt.maxrandomp) nrand=1
555          if (nrand+ifine.gt.maxrand2) nrand=1
556          rhoaux=rhograd/rhoa
557          dtf=dt*fine
558
559          dtftlw=dtf/tlw
560
561! Loop over ifine short time steps for vertical component
562!********************************************************
563
564          do i=1,ifine
565
566! Determine the drift velocity and density correction velocity
567! Determine the drift velocity and density correction velocity
568!*************************************************************
569!--------------- lines below are teh original FLEXPART  and are commented out to insert teh cbl options comment by mc
570!            if (turbswitch) then
571!              if (dtftlw.lt..5) then
572!                wp=((1.-dtftlw)*wp+rannumb(nrand+i)*sqrt(2.*dtftlw) &
573!                +dtf*(dsigwdz+rhoaux*sigw))*real(icbt)
574!              else
575!                rw=exp(-dtftlw)
576!                wp=(rw*wp+rannumb(nrand+i)*sqrt(1.-rw**2) &
577!                +tlw*(1.-rw)*(dsigwdz+rhoaux*sigw))*real(icbt)
578!              endif
579!              delz=wp*sigw*dtf
580!            else
581!              rw=exp(-dtftlw)
582!              wp=(rw*wp+rannumb(nrand+i)*sqrt(1.-rw**2)*sigw &
583!              +tlw*(1.-rw)*(dsigw2dz+rhoaux*sigw**2))*real(icbt)
584!              delz=wp*dtf
585!            endif
586!************ CBL scheme integrated in FLEXPART: added by mc **********!
587         if (turbswitch) then
588          if (dtftlw.lt..5) then
589            if (cblflag.eq.1) then
590!     wp2=wp
591!     zt2=zt
592!    ust2=ust
593!    wst2=wst
594!    h2=h
595!    rhoa2=rhoa
596!    rhograd2=rhograd
597!    sigw2=sigw
598!    dsigwdz2=dsigwdz
599!    tlw2=tlw
600!    ptot_lhh2=ptot_lhh
601!    Q_lhh2=Q_lhh
602!    phi_lhh2=phi_lhh
603!    ath2=ath
604!    bth2=bth
605!    ol2=ol
606               if (-h/ol.gt.5) then  !modified by mc
607              !if (ol.lt.0.) then   !modified by mc 
608              !if (ol.gt.0.) then   !modified by mc : for test
609                  !print  *,zt,wp,ath,bth,tlw,dtf,'prima'
610                  flagrein=0
611                  nrand=nrand+1
612                  old_wp_buf=wp
613                  del_test=(1.-old_wp_buf)/old_wp_buf
614
615                  !rhoa=1.  !for testing vertical well mixed state, by mc
616                  !rhograd=0.  !for testing vertical well mixed state, by mc
617                  call cbl(wp,zt,ust,wst,h,rhoa,rhograd,sigw,dsigwdz,tlw,ptot_lhh,Q_lhh,phi_lhh,ath,bth,ol,flagrein)
618! see inside the routine for inverse time
619                  wp=(wp+ath*dtf+bth*rannumb(nrand)*sqrt(dtf))*real(icbt)
620                  delz=wp*dtf
621!                  if ((ieee_is_nan(zt).or.ieee_is_nan(wp)).and.(flagrein.eq.0)) print*,'pb4',wp2,zt2,ust2,wst2,h2,rhoa2,rhograd2,sigw2,dsigwdz2,tlw2,ptot_lhh2,Q_lhh2,phi_lhh2,ath2,bth2,ol2,flagrein,i
622                  if (abs(wp).gt.50.) flagrein=1
623                  if (flagrein.eq.1) then  !added for re-initlization of particle vertical velcoity based on condition inside routine cbl.f90
624                      call re_initialize_particle(zt,ust,wst,h,sigw,old_wp_buf,nrand,ol)
625!                  if (ieee_is_nan(old_wp_buf)) print*,"PROBLEM WP",wp,old_wp_buf,nrand,ol,zt,ust,wst,h,sigw
626                      wp=old_wp_buf
627                      delz=wp*dtf
628                      !nan_count(myid)=nan_count(myid)+1
629                       nan_count(ompid+1)=nan_count(ompid+1)+1
630                  else       
631                  del_test=(1.-wp)/wp !catch infinity value
632!                 if (ieee_is_nan(wp) .or. ieee_is_nan(del_test)) then
633                   if (isnan2(wp).or.isnan2(del_test)) then !note that, given the test on particle velocity inside the routine cbl.f90, this condition should never be true!!
634!                 if (isnan(wp).or.isnan(del_test)) then !note that, given the test on particle velocity inside the routine cbl.f90, this condition should never be true!!
635                      nrand=nrand+1
636                      call re_initialize_particle(zt,ust,wst,h,sigw,old_wp_buf,nrand,ol)
637                      wp=old_wp_buf
638                      delz=wp*dtf
639                      nan_count(ompid+1)=nan_count(ompid+1)+1
640                      !nan_count(myid)=nan_count(myid)+1
641                      print *,'NaN counter equal to:',nan_count(ompid+1),'omp',ompid,'mpi',myid
642!                     ,'increase
643!                     !ifine if this number became a non-negligible fraction of
644!                     !the particle number'
645                  end if
646                  end if
647              else
648                  !rhoa=1.  !for testing vertical well mixed state, by mc
649                  !rhograd=0. !for testing vertical well mixed state, by mc
650                  nrand=nrand+1
651                  old_wp_buf=wp                 
652                  ath=-wp/tlw+sigw*dsigwdz+wp*wp/sigw*dsigwdz+sigw*sigw/rhoa*rhograd  !1-note for inverse time should be -wp/tlw*ldirect+... calculated for wp=-wp
653                                                                                      !2-but since ldirect =-1 for inverse time and this must be calculated for (-wp) and
654                                                                                      !3-the gaussian pdf is symmetric (i.e. pdf(w)=pdf(-w) ldirect can be discarded
655                  bth=sigw*rannumb(nrand)*sqrt(2.*dtftlw)
656                  wp=(wp+ath*dtf+bth)*real(icbt) 
657                  delz=wp*dtf
658                  del_test=(1.-wp)/wp !catch infinity value
659!                    if (ieee_is_nan(wp).or.ieee_is_nan(del_test)) then
660!                print*,'PB',wp2,zt2,ust2,wst2,h2,rhoa2,rhograd2,sigw2,dsigwdz2,tlw2,ptot_lhh2,Q_lhh2,phi_lhh2,ath2,bth2,ol2,flagrein,i
661!                print*,'PB2',ath,old_wp_buf,bth,wp,sigw
662                  if (isnan2(wp).or.isnan2(del_test).or.abs(wp).gt.50.) then
663!                  if (wp.ne.wp .or. del_test.ne.del_test) then
664!                  if (ieee_is_nan(wp) .or. ieee_is_nan(del_test).or.abs(wp).gt.50.) then           
665                      nrand=nrand+1                     
666                      wp=sigw*rannumb(nrand)
667                      delz=wp*dtf
668                      nan_count(ompid+1)=nan_count(ompid+1)+1
669                      print *,'NaN counter equal to:',nan_count(ompid+1),'omp',ompid,'mpi',myid &
670            ,'increase ifine if this number became a non-negligible fraction of the particle number'
671                  end if 
672              end if
673            else
674                 wp=((1.-dtftlw)*wp+rannumb(nrand+i)*sqrt(2.*dtftlw) &
675                 +dtf*(dsigwdz+rhoaux*sigw))*real(icbt)
676                 delz=wp*sigw*dtf
677            end if
678          else
679            rw=exp(-dtftlw)
680            wp=(rw*wp+rannumb(nrand+i)*sqrt(1.-rw**2) &
681                 +tlw*(1.-rw)*(dsigwdz+rhoaux*sigw))*real(icbt)
682            delz=wp*sigw*dtf
683          endif
684         
685        else
686          rw=exp(-dtftlw)
687          wp=(rw*wp+rannumb(nrand+i)*sqrt(1.-rw**2)*sigw &
688               +tlw*(1.-rw)*(dsigw2dz+rhoaux*sigw**2))*real(icbt)
689          delz=wp*dtf
690        endif
691
692!***************** end turbulent options : comemnt by mc *********************************!
693
694!        if (ieee_is_nan(wp)) then
695!         print*,"PROBLEM WP OUT",wp,old_wp_buf,nrand,ol,zt,ust,wst,h,sigw
696!        endif
697! FLEXPART_WRF - zero up,vp,wp if turbulence is turned off
698            if (turb_option .eq. turb_option_none) then
699             up=0.0
700             vp=0.0
701              wp=0.0
702              delz=0.
703            end if
704!          print*,'delz',delz,zt
705!****************************************************
706! Compute turbulent vertical displacement of particle
707!****************************************************
708!     if (nrelpoint.eq.970) then
709!      write(15,*),rannumb(nrand),nrand,nrelpoint,OMP_GET_THREAD_NUM()
710!     endif
711
712            if (abs(delz).gt.h) delz=mod(delz,h)
713
714! Determine if particle transfers to a "forbidden state" below the ground
715! or above the mixing height
716!************************************************************************
717
718            if (delz.lt.-zt) then         ! reflection at ground
719              icbt=-1
720              zt=-zt-delz
721               else if (delz.gt.(h-zt)) then ! reflection at h
722                 icbt=-1
723                 zt=-zt-delz+2.*h
724!            else if (delz.gt.(h-zt) .and. reflect_switch==1) then ! reflection at h
725!             else if (delz.gt.(h-zt)) then ! reflection at h
726!               icbt=-1
727!               zt=-zt-delz+2.*h
728            else                         ! no reflection
729              icbt=1
730              zt=zt+delz
731            endif
732
733            if (i.ne.ifine) then
734! FLEXPART_WRF, TKE option
735              if (turb_option .gt. 1) then
736                  do ii=2,nz
737                      if (height(ii).gt.zt) then
738                       indz=ii-1
739                       indzp=ii
740                       goto 69
741                      endif
742                   enddo
74369                 continue
744 
745! If one of the levels necessary is not yet available,
746! calculate it
747!*****************************************************
748 
749            do ii=indz,indzp                !i
750              if (indzindicator(ii)) then
751                if (ngrid.le.0) then
752                  call interpol_misslev(ii,xt,yt,zt, &
753    uprof,vprof,wprof, usigprof,vsigprof,wsigprof, &
754    rhoprof,rhogradprof, tkeprof,pttprof, &
755    u,v,w,usig,vsig,wsig,pvi, &
756    p1,p2,p3,p4,ddx,ddy,rddx,rddy,dtt,dt1,dt2, &
757    ix,jy,ixp,jyp,ngrid,indz,indzp, depoindicator, &
758    indzindicator, &
759    ust,wst,ol,h,zeta,sigu,sigv,tlu,tlv,tlw, &
760    sigw,dsigwdz,dsigw2dz)
761                else
762                  call interpol_misslev_nests(ii,xt,yt,zt, &
763    uprof,vprof,wprof, usigprof,vsigprof,wsigprof, &
764    rhoprof,rhogradprof, tkeprof,pttprof, &
765    u,v,w,usig,vsig,wsig,pvi, &
766    p1,p2,p3,p4,ddx,ddy,rddx,rddy,dtt,dt1,dt2, &
767    ix,jy,ixp,jyp,ngrid,indz,indzp, depoindicator, &
768    indzindicator, &
769    ust,wst,ol,h,zeta,sigu,sigv,tlu,tlv,tlw, &
770    sigw,dsigwdz,dsigw2dz)
771                endif
772              endif
773            enddo                          !i
774!              write(*,*)'after reflection'
775            if(turb_option .eq. turb_option_mytke) &
776               call tke_partition_my(zt, &
777   ust,wst,ol,h,zeta,sigu,sigv,tlu,tlv,tlw, &
778   sigw,dsigwdz,dsigw2dz,uprof,vprof,tkeprof,pttprof,indz,indzp)
779            if(turb_option .eq. turb_option_tke)   &
780               call tke_partition_hanna(zt, &
781   ust,wst,ol,h,zeta,sigu,sigv,tlu,tlv,tlw, &
782   sigw,dsigwdz,dsigw2dz,uprof,vprof,tkeprof,pttprof,indz,indzp)
783            else
784                zeta=zt/h
785                 call hanna_short(zt, &
786    ust,wst,ol,h,zeta,sigu,sigv,tlu,tlv,tlw, &
787    sigw,dsigwdz,dsigw2dz)
788
789            endif
790           endif
791
792          enddo
793         
794          if (cblflag.ne.1) nrand=nrand+i  !------>>>>>>>>>>>>>>>> modified by mc for accounting of different increment of nrand in cbl flag
795!    if (nombre.eq.103) print*,'usig 3',usig
796
797! Determine time step for next integration
798!*****************************************
799
800          if (turbswitch) then
801            ldt=int(min(tlw,h/max(2.*abs(wp*sigw),1.e-5), &
802            0.5/abs(dsigwdz))*ctl)
803          else
804            ldt=int(min(tlw,h/max(2.*abs(wp),1.e-5))*ctl)
805          endif
806          ldt=max(ldt,mintime)
807
808
809! If particle represents only a single species, add gravitational settling
810! velocity. The settling velocity is zero for gases, or if particle
811! represents more than one species
812!*************************************************************************
813
814      if (mdomainfill.eq.0) then
815        do nsp=1,nspec
816!         print*,nrelpoint,nsp
817          if (xmass(nrelpoint,nsp).gt.eps2) goto 887
818        end do
819887     nsp=min(nsp,nspec)
820       if (density(nsp).gt.0.) then
821!        print*,'settle'
822!       print*,'settle 1'
823               call get_settling(itime,real(xt),real(yt),zt,nsp,settling)  !bugfix
824        endif
825        w=w+settling
826      endif
827
828
829! Horizontal displacements during time step dt are small real values compared
830! to the position; adding the two, would result in large numerical errors.
831! Thus, displacements are accumulated during lsynctime and are added to the
832! position at the end
833!****************************************************************************
834
835          dxsave=dxsave+u*dt
836!     if (nombre.eq.103) print*,'xt-2',dxsave,u,dt
837          dysave=dysave+v*dt
838          dawsave=dawsave+up*dt
839          dcwsave=dcwsave+vp*dt
840          zt=zt+w*dt*real(ldirect)  ! comment out and put zt=zt for testing equation based on the well_mixed conditin comemnt by mc
841
842
843          if (zt.gt.h) then
844            if (itimec.eq.itime+lsynctime) goto 99
845            goto 700    ! complete the current interval above PBL
846          endif
847          if (zt.lt.0.) zt=-1.*zt    ! if particle below ground -> refletion
848
849
850!!!! CHANGE: TEST OF THE WELL-MIXED CRITERION
851!!!! These lines may be switched on to test the well-mixed criterion
852!      if (zt.le.h) then
853!        zacc=zacc+zt/h*dt
854!        hsave=hsave+h*dt
855!        tacc=tacc+dt
856!        do 67 i=1,iclass
857!          if ((zt/h.gt.th(i-1)).and.(zt/h.le.th(i)))
858!     +    t(i)=t(i)+dt
859!67        continue
860!      endif
861!c       print*,'itime',itime
862!c     if ((mod(abs(itime),3600).eq.0)) then
863!c     if ((mod(abs(itime),3600).eq.0).and.dump) then
864!      if (itime<itimeold) then
865!       print*,'dump well',itime,itimeold
866!       dump=.false.
867!        itimeold=itimeold-3600
868!       write(550,'(i8,12f10.3)') itime,hsave/tacc,zacc/tacc,
869!c      write(550,'(i8,22f10.3)') itime,hsave/tacc,zacc/tacc,
870!     + (t(i)/tacc*real(iclass),i=1,iclass)
871!        flush(550)
872!        zacc=0.
873!        tacc=0.
874!        do 68 i=1,iclass
875!68        t(i)=0.
876!        hsave=0.
877!      endif
878!      if (mod(abs(itime),3600).ne.0) dump=.true.
879!c       print*,'itime',itime,3600,mod(abs(itime),3600),dump
880!!!! CHANGE
881!!!******************  NEW test for THE WELL MIXED CRITERION by mc ***************
882!$OMP CRITICAL
883!if (zt.lt.h) then
884!          i_well=int(zt/h*num_lvl*1.)+1                                          !per fare il test qui devo considerare OMP and MPI...
885!          well_mixed_vector(i_well,ompid+1)=well_mixed_vector(i_well,ompid+1)+dt
886!          well_mixed_norm(ompid+1)=well_mixed_norm(ompid+1)+dt
887!          avg_air_dens(i_well,ompid+1)=avg_air_dens(i_well,ompid+1)+rhoa*dt
888!         
889!      end if
890!      h_well(ompid+1)=h
891!$OMP END CRITICAL   
892!!*********************************************************************************
893
894! Determine probability of deposition
895!************************************
896
897      if ((DRYDEP).and.(zt.lt.2.*href)) then
898        do ks=1,nspec
899          if (DRYDEPSPEC(ks)) then
900            if (depoindicator(ks)) then
901              if (ngrid.le.0) then
902                   call interpol_vdep(ks,vdepo(ks),ix,jy,ixp,jyp, &
903                                      p1,p2,p3,p4,dt1,dt2,dtt,depoindicator)
904              else
905                call interpol_vdep_nests(ks,vdepo(ks),ix,jy,ixp,jyp, &
906                                      p1,p2,p3,p4,dt1,dt2,dtt,depoindicator,ngrid)
907              endif
908            endif
909  ! correction by Petra Seibert, 10 April 2001
910  !   this formulation means that prob(n) = 1 - f(0)*...*f(n)
911  !   where f(n) is the exponential term
912               prob(ks)=1.+(prob(ks)-1.)* &
913                    exp(-vdepo(ks)*abs(dt)/(2.*href))
914          endif
915        end do
916      endif
917
918      if (zt.lt.0.) zt=min(h-eps2,-1.*zt)    ! if particle below ground -> reflection
919
920      if (itimec.eq.(itime+lsynctime)) then
921!   if (nombre.eq.103) print*,'usig',usig,usigprof(indzp)+usigprof(indz),indz
922        usig=0.5*(usigprof(indzp)+usigprof(indz))
923        vsig=0.5*(vsigprof(indzp)+vsigprof(indz))
924        wsig=0.5*(wsigprof(indzp)+wsigprof(indz))
925        goto 99  ! finished
926      endif
927      goto 100
928
929! END TIME LOOP
930!==============
931
932
933      endif
934
935
936
937!**********************************************************
938! For all particles that are outside the PBL, make a single
939! time step. Only horizontal turbulent disturbances are
940! calculated. Vertical disturbances are reset.
941!**********************************************************
942
943
944! Interpolate the wind
945!*********************
946
947!JB needs to define mu and mv
948700   continue
949      if (ngrid.le.0) then
950        xts=real(xt)
951        yts=real(yt)
952        call interpol_wind(itime,xts,yts,zt, &
953    uprof,vprof,wprof, usigprof,vsigprof,wsigprof, &
954    rhoprof,rhogradprof, tkeprof,pttprof, &
955    u,v,w,usig,vsig,wsig,pvi, &
956    p1,p2,p3,p4,ddx,ddy,rddx,rddy,dtt,dt1,dt2, &
957    ix,jy,ixp,jyp,ngrid,indz,indzp, depoindicator, &
958    indzindicator,mu,mv)
959 !JB mw not needed here
960      else
961        call interpol_wind_nests(itime,xtn,ytn,zt, &
962    uprof,vprof,wprof, usigprof,vsigprof,wsigprof, &
963    rhoprof,rhogradprof, tkeprof,pttprof, &
964    u,v,w,usig,vsig,wsig,pvi, &
965    p1,p2,p3,p4,ddx,ddy,rddx,rddy,dtt,dt1,dt2, &
966    ix,jy,ixp,jyp,ngrid,indz,indzp, depoindicator, &
967    indzindicator,mu,mv)
968      endif
969
970!    if (nombre.eq.103) print*,'usig 4',usig
971! Compute everything for above the PBL
972
973! Assume constant, uncorrelated, turbulent perturbations
974! In the stratosphere, use a small vertical diffusivity d_strat,
975! in the troposphere, use a larger horizontal diffusivity d_trop.
976! Turbulent velocity scales are determined based on sqrt(d_trop/dt)
977!******************************************************************
978
979      ldt=abs(lsynctime-itimec+itime)
980      dt=real(ldt)
981     
982  if (zt.lt.tropop) then  ! in the troposphere
983    uxscale=sqrt(2.*d_trop/dt)
984!   if (nrand+1.gt.maxrandomp) nrand=1
985    if (nrand+1.gt.maxrand2) nrand=1
986    ux=rannumb(nrand)*uxscale
987    vy=rannumb(nrand+1)*uxscale
988    nrand=nrand+2
989    wp=0.
990  else if (zt.lt.tropop+1000.) then     ! just above the tropopause: make transition
991    weight=(zt-tropop)/1000.
992    uxscale=sqrt(2.*d_trop/dt*(1.-weight))
993!   if (nrand+2.gt.maxrandomp) nrand=1
994    if (nrand+2.gt.maxrand2) nrand=1
995    ux=rannumb(nrand)*uxscale
996    vy=rannumb(nrand+1)*uxscale
997    wpscale=sqrt(2.*d_strat/dt*weight)
998    wp=rannumb(nrand+2)*wpscale+d_strat/1000.
999    nrand=nrand+3
1000  else                 ! in the stratosphere
1001!   if (nrand.gt.maxrandomp) nrand=1
1002    if (nrand.gt.maxrand2) nrand=1
1003    ux=0.
1004    vy=0.
1005    wpscale=sqrt(2.*d_strat/dt)
1006    wp=rannumb(nrand)*wpscale
1007    nrand=nrand+1
1008  endif
1009
1010! FLEXPART_WRF - zero ux,vy,wp if turbulence is turned off
1011      if (turb_option .eq. turb_option_none) then
1012        ux=0.0
1013        vy=0.0
1014        wp=0.0
1015      end if
1016
1017
1018! If particle represents only a single species, add gravitational settling
1019! velocity. The settling velocity is zero for gases
1020!*************************************************************************
1021
1022    if (mdomainfill.eq.0) then
1023      do nsp=1,nspec
1024        if (xmass(nrelpoint,nsp).gt.eps2) goto 888
1025      end do
1026888   nsp=min(nsp,nspec)
1027    if (density(nsp).gt.0.) then
1028!       print*,'settle 2, bef',real(xt),real(yt),zt,cpt
1029           call get_settling(itime,real(xt),real(yt),zt,nsp,settling)  !bugfix
1030!       print*,'settle 2, aft',real(xt),real(yt),zt,cpt,w
1031!        print*,'settle'
1032      endif
1033      w=w+settling
1034    endif
1035
1036! Calculate position at time step itime+lsynctime
1037!************************************************
1038
1039!       print*,'settle 2, aft1.5',zt,settling,wp,dt
1040      dxsave=dxsave+(u+ux)*dt
1041!     if (nombre.eq.103) print*,'xt-1',dxsave,u,ux,dt
1042      dysave=dysave+(v+vy)*dt
1043      zt=zt+(w+wp)*dt*real(ldirect)
1044!       print*,'settle 2, aft2',zt,cpt
1045      if (zt.lt.0.) zt=min(h-eps2,-1.*zt)    ! if particle below ground -> reflection
1046!       print*,'settle 2, aft3',zt,cpt
1047
104899    continue
1049
1050
1051
1052!****************************************************************
1053! Add mesoscale random disturbances
1054! This is done only once for the whole lsynctime interval to save
1055! computation time
1056!****************************************************************
1057
1058
1059! Mesoscale wind velocity fluctuations are obtained by scaling
1060! with the standard deviation of the grid-scale winds surrounding
1061! the particle location, multiplied by a factor turbmesoscale.
1062! The autocorrelation time constant is taken as half the
1063! time interval between wind fields
1064!****************************************************************
1065
1066      r=exp(-2.*real(abs(lsynctime))/real(lwindinterv))
1067      rs=sqrt(1.-r**2)
1068!     if (nrand+2.gt.maxrandomp) nrand=1
1069      if (nrand+2.gt.maxrand2) nrand=1
1070!     if (nombre.eq.103) print*,'usgig0',r,usigold,rannumb(nrand)
1071      usigold=r*usigold+rs*rannumb(nrand)*usig*turbmesoscale
1072!     if (nombre.eq.103) print*,'usgig1',usigold,usig,turbmesoscale
1073      vsigold=r*vsigold+rs*rannumb(nrand+1)*vsig*turbmesoscale
1074      wsigold=r*wsigold+rs*rannumb(nrand+2)*wsig*turbmesoscale
1075
1076! FLEXPART_WRF - zero u/v/wsigold if turbulence is turned off
1077! Note: for mesoscale model applications this component should be ignored!
1078!     if (turb_option .eq. turb_option_none) then
1079!       usigold=0.0
1080!       vsigold=0.0
1081!       wsigold=0.0
1082!     end if
1083
1084      dxsave=dxsave+usigold*real(lsynctime)
1085      dysave=dysave+vsigold*real(lsynctime)
1086
1087      zt=zt+wsigold*real(lsynctime)
1088!       print*,'settle 2, aft4',zt,cpt
1089      if (zt.lt.0.) zt=-1.*zt    ! if particle below ground -> refletion
1090!       print*,'settle 2, aft5',zt,cpt
1091
1092!*************************************************************
1093! Transform along and cross wind components to xy coordinates,
1094! add them to u and v, transform u,v to grid units/second
1095! and calculate new position
1096!*************************************************************
1097
1098      call windalign(dxsave,dysave,dawsave,dcwsave,ux,vy)
1099      dxsave=dxsave+ux
1100!     if (nombre.eq.103) print*,'xt0',dxsave,usigold,ux
1101      dysave=dysave+vy
1102      if (ngrid.ge.0) then
1103! for FLEXPART_WRF, dx & dy are in meters,
1104! dxconst=1/dx, dyconst=1/dy, and no cos(lat) is needed
1105!       cosfact=dxconst/cos((yt*dy+ylat0)*pi180)
1106!     if (nombre.eq.103) print*,'xt1',xt,dxsave,dxconst
1107!       xt=xt+real(dxsave*dxconst*real(ldirect),kind=dp)
1108!       yt=yt+real(dysave*dyconst*real(ldirect),kind=dp)
1109!      xt=xt+real(dxsave/mu*dxconst*real(ldirect),kind=dp)
1110!      yt=yt+real(dysave/mv*dyconst*real(ldirect),kind=dp)
1111       xt=xt +real(dxsave*mu*dxconst*real(ldirect),kind=dp)  !IF COOMMENTED OUT TO is to ISOLate VERTCAL FORMULAITON FOR TEST REASON BY mc
1112       yt=yt +real(dysave*mv*dyconst*real(ldirect),kind=dp)  !IF COOMMENTED OUT TO is to ISOLate VERTCAL FORMULAITON FOR TEST REASON BY mc
1113! JB: needs interpolate m_w on the coordinates
1114!     else if (ngrid.eq.-1) then      ! around north pole
1115!       xlon=xlon0+xt*dx
1116!       ylat=ylat0+yt*dy
1117!       call cll2xy(northpolemap,ylat,xlon,xpol,ypol)
1118!       gridsize=1000.*cgszll(northpolemap,ylat,xlon)
1119!       dxsave=dxsave/gridsize
1120!       dysave=dysave/gridsize
1121!       xpol=xpol+dxsave*real(ldirect)
1122!       ypol=ypol+dysave*real(ldirect)
1123!       call cxy2ll(northpolemap,xpol,ypol,ylat,xlon)
1124!       xt=(xlon-xlon0)/dx
1125!       yt=(ylat-ylat0)/dy
1126!     else if (ngrid.eq.-2) then    ! around south pole
1127!       xlon=xlon0+xt*dx
1128!       ylat=ylat0+yt*dy
1129!       call cll2xy(southpolemap,ylat,xlon,xpol,ypol)
1130!       gridsize=1000.*cgszll(southpolemap,ylat,xlon)
1131!       dxsave=dxsave/gridsize
1132!       dysave=dysave/gridsize
1133!       xpol=xpol+dxsave*real(ldirect)
1134!       ypol=ypol+dysave*real(ldirect)
1135!       call cxy2ll(southpolemap,xpol,ypol,ylat,xlon)
1136!       xt=(xlon-xlon0)/dx
1137!       yt=(ylat-ylat0)/dy
1138      else
1139        write(*,*) 'advance -- bad ngrid = ', ngrid
1140        stop
1141      endif
1142
1143
1144! If global data are available, use cyclic boundary condition
1145!************************************************************
1146
1147      if (xglobal) then
1148        if (xt.ge.real(nxmin1)) xt=xt-real(nxmin1)
1149        if (xt.lt.0.) xt=xt+real(nxmin1)
1150        if (xt.le.eps) xt=eps
1151        if (abs(xt-real(nxmin1)).le.eps) xt=real(nxmin1)-eps
1152      endif
1153
1154
1155! Check position: If trajectory outside model domain, terminate it
1156!*****************************************************************
1157
1158      if ((xt.lt.0.).or.(xt.ge.real(nxmin1)).or.(yt.lt.0.).or. &
1159      (yt.ge.real(nymin1))) then
1160        nstop=3
1161        return
1162      endif
1163
1164! If particle above highest model level, set it back into the domain
1165!*******************************************************************
1166
1167      if (zt.ge.height(nz)) zt=height(nz)-100.*eps
1168!       print*,'settle 2, aft6',zt,cpt
1169
1170
1171!************************************************************************
1172! Now we could finish, as this was done in FLEXPART versions up to 4.0.
1173! However, truncation errors of the advection can be significantly
1174! reduced by doing one iteration of the Petterssen scheme, if this is
1175! possible.
1176! Note that this is applied only to the grid-scale winds, not to
1177! the turbulent winds.
1178!************************************************************************
1179
1180! The Petterssen scheme can only applied with long time steps (only then u
1181! is the "old" wind as required by the scheme); otherwise do nothing
1182!*************************************************************************
1183
1184      if (ldt.ne.abs(lsynctime)) return
1185
1186! The Petterssen scheme can only be applied if the ending time of the time step
1187! (itime+ldt*ldirect) is still between the two wind fields held in memory;
1188! otherwise do nothing
1189!******************************************************************************
1190
1191      if (abs(itime+ldt*ldirect).gt.abs(memtime(2))) return 
1192
1193! Apply it also only if starting and ending point of current time step are on
1194! the same grid; otherwise do nothing
1195!*****************************************************************************
1196      if (nglobal.and.(yt.gt.switchnorthg)) then
1197        ngr=-1
1198      else if (sglobal.and.(yt.lt.switchsouthg)) then
1199        ngr=-2
1200      else
1201    ngr=0
1202    do j=numbnests,1,-1
1203      if ((xt.gt.xln(j)+eps).and.(xt.lt.xrn(j)-eps).and. &
1204           (yt.gt.yln(j)+eps).and.(yt.lt.yrn(j)-eps)) then
1205        ngr=j
1206        goto 43
1207      endif
1208    end do
120943   continue
1210  endif
1211
1212      if (ngr.ne.ngrid) return
1213
1214! Determine nested grid coordinates
1215!**********************************
1216
1217      if (ngrid.gt.0) then
1218        xtn=(xt-xln(ngrid))*xresoln(ngrid)
1219        ytn=(yt-yln(ngrid))*yresoln(ngrid)
1220        ix=int(xtn)
1221        jy=int(ytn)
1222      else
1223        ix=int(xt)
1224        jy=int(yt)
1225      endif
1226      ixp=ix+1
1227      jyp=jy+1
1228
1229
1230! Memorize the old wind
1231!**********************
1232
1233      uold=u
1234      vold=v
1235      wold=w
1236
1237! Interpolate wind at new position and time
1238!******************************************
1239
1240      if (ngrid.le.0) then
1241        xts=real(xt)
1242        yts=real(yt)
1243        call interpol_wind_short(itime+ldt*ldirect,xts,yts,zt, &
1244    uprof,vprof,wprof, usigprof,vsigprof,wsigprof, &
1245    rhoprof,rhogradprof, tkeprof,pttprof, &
1246    u,v,w,usig,vsig,wsig,pvi, &
1247    p1,p2,p3,p4,ddx,ddy,rddx,rddy,dtt,dt1,dt2, &
1248    ix,jy,ixp,jyp,ngrid,indz,indzp, depoindicator, &
1249    indzindicator)
1250      else
1251        call interpol_wind_short_nests(itime+ldt*ldirect,xtn,ytn,zt, &
1252    uprof,vprof,wprof, usigprof,vsigprof,wsigprof, &
1253    rhoprof,rhogradprof, tkeprof,pttprof, &
1254    u,v,w,usig,vsig,wsig,pvi, &
1255    p1,p2,p3,p4,ddx,ddy,rddx,rddy,dtt,dt1,dt2, &
1256    ix,jy,ixp,jyp,ngrid,indz,indzp, depoindicator, &
1257    indzindicator)
1258
1259      endif
1260!       print*,'settle 2, aft7',zt,cpt
1261
1262  if (mdomainfill.eq.0) then
1263    do nsp=1,nspec
1264      if (xmass(nrelpoint,nsp).gt.eps2) goto 889
1265    end do
1266889   nsp=min(nsp,nspec)
1267   if (density(nsp).gt.0.) then
1268!       print*,'settle 3, bef',real(xt),real(yt),zt,cpt
1269        call get_settling(itime+ldt,real(xt),real(yt),zt,nsp,settling)  !bugfix
1270!       print*,'settle 3, aft',real(xt),real(yt),zt,cpt
1271!        print*,'settle'
1272      endif
1273    w=w+settling
1274  endif
1275
1276! Determine the difference vector between new and old wind
1277! (use half of it to correct position according to Petterssen)
1278!*************************************************************
1279
1280      u=(u-uold)/2.
1281      v=(v-vold)/2.
1282      w=(w-wold)/2.
1283
1284
1285! Finally, correct the old position
1286!**********************************
1287
1288      zt=zt+w*real(ldt*ldirect)
1289  if (zt.lt.0.) zt=min(h-eps2,-1.*zt)    ! if particle below ground -> reflection
1290
1291      if (ngrid.ge.0) then
1292! for FLEXPART_WRF, dx & dy are in meters,
1293! dxconst=1/dx, dyconst=1/dy, and no cos(lat) is needed
1294!       cosfact=dxconst/cos((yt*dy+ylat0)*pi180)
1295!       if (nombre.eq.103)  print*,'xt2',xt,u,dxconst,ldt
1296!       xt=xt+real(u*dxconst*real(ldt*ldirect),kind=dp)
1297!       yt=yt+real(v*dyconst*real(ldt*ldirect),kind=dp)
1298!        print*,'mw',mu,mv
1299!        xt=xt+real(u*dxconst/mu*real(ldt*ldirect),kind=dp)
1300!        yt=yt+real(v*dyconst/mv*real(ldt*ldirect),kind=dp)
1301         xt=xt +real(u*dxconst*mu*real(ldt*ldirect),kind=dp)  !IF COOMMENTED OUT TO is to ISOLate VERTCAL FORMULAITON FOR TEST REASON BY mc
1302         yt=yt +real(v*dyconst*mv*real(ldt*ldirect),kind=dp)  !IF COOMMENTED OUT TO is to ISOLate VERTCAL FORMULAITON FOR TEST REASON BY mc
1303!     else if (ngrid.eq.-1) then      ! around north pole
1304!       xlon=xlon0+xt*dx
1305!       ylat=ylat0+yt*dy
1306!       call cll2xy(northpolemap,ylat,xlon,xpol,ypol)
1307!       gridsize=1000.*cgszll(northpolemap,ylat,xlon)
1308!       u=u/gridsize
1309!       v=v/gridsize
1310!       xpol=xpol+u*real(ldt*ldirect)
1311!       ypol=ypol+v*real(ldt*ldirect)
1312!       call cxy2ll(northpolemap,xpol,ypol,ylat,xlon)
1313!       xt=(xlon-xlon0)/dx
1314!       yt=(ylat-ylat0)/dy
1315!     else if (ngrid.eq.-2) then    ! around south pole
1316!       xlon=xlon0+xt*dx
1317!       ylat=ylat0+yt*dy
1318!       call cll2xy(southpolemap,ylat,xlon,xpol,ypol)
1319!       gridsize=1000.*cgszll(southpolemap,ylat,xlon)
1320!       u=u/gridsize
1321!       v=v/gridsize
1322!       xpol=xpol+u*real(ldt*ldirect)
1323!       ypol=ypol+v*real(ldt*ldirect)
1324!       call cxy2ll(southpolemap,xpol,ypol,ylat,xlon)
1325!       xt=(xlon-xlon0)/dx
1326!       yt=(ylat-ylat0)/dy
1327      else
1328        write(*,*) 'advance -- bad ngrid = ', ngrid
1329        stop
1330      endif
1331
1332! If global data are available, use cyclic boundary condition
1333!************************************************************
1334
1335      if (xglobal) then
1336        if (xt.ge.real(nxmin1)) xt=xt-real(nxmin1)
1337        if (xt.lt.0.) xt=xt+real(nxmin1)
1338        if (xt.le.eps) xt=eps
1339        if (abs(xt-real(nxmin1)).le.eps) xt=real(nxmin1)-eps
1340      endif
1341
1342! Check position: If trajectory outside model domain, terminate it
1343!*****************************************************************
1344
1345      if ((xt.lt.0.).or.(xt.ge.real(nxmin1)).or.(yt.lt.0.).or. &
1346      (yt.ge.real(nymin1))) then
1347        nstop=3
1348        return
1349      endif
1350
1351! If particle above highest model level, set it back into the domain
1352!*******************************************************************
1353
1354      if (zt.ge.height(nz)) zt=height(nz)-100.*eps
1355
1356!       if (nombre.eq.103)  print*,'end',xt,u,dxconst,ldt
1357
1358      end subroutine advance
1359
1360!       logical function isnan2(a)
1361!       real :: a
1362!        if ((a.ne.a)) then !.or.((a*0.).ne.0.)) then
1363!       isnan2 = .true.
1364!       else
1365!       isnan2 = .false.
1366!        end if
1367!        return
1368!       end
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG