Changeset 20 for trunk/src/readreleases.f90

Timestamp:

Dec 23, 2013, 6:23:38 PM (11 years ago)

Author:

igpis

Message:

move version 9.1.8 form branches to trunk. Contributions from HSO, saeck, pesei, NIK, RT, XKF, IP and others

File:

• trunk/src/readreleases.f90 (modified) (5 diffs)

trunk/src/readreleases.f90

@@ -57 +57 @@
   ! weta, wetb          parameters to determine the wet scavenging coefficient *
   ! zpoint1,zpoint2     height range, over which release takes place           *
+  ! num_min_discrete    if less, release cannot be randomized and happens at   *
+  !                     time mid-point of release interval                     *
   !                                                                            *
   !*****************************************************************************
@@ -68 +70 @@
 
   integer :: numpartmax,i,j,id1,it1,id2,it2,specnum_rel,idum,stat
+  integer,parameter :: num_min_discrete=100
   real :: vsh(ni),fracth(ni),schmih(ni),releaserate,xdum,cun
-  real(kind=dp) :: jul1,jul2,juldate
+  real(kind=dp) :: jul1,jul2,julm,juldate
   character(len=50) :: line
   logical :: old
@@ -376 +379 @@
   jul1=juldate(id1,it1)
   jul2=juldate(id2,it2)
+  julm=(jul1+jul2)/2.
   if (jul1.gt.jul2) then
     write(*,*) 'FLEXPART MODEL ERROR'
@@ -391 +395 @@
         stop
       endif
-      ireleasestart(numpoint)=int((jul1-bdate)*86400.)
-      ireleaseend(numpoint)=int((jul2-bdate)*86400.)
+      if (npart(numpoint).gt.num_min_discrete) then
+        ireleasestart(numpoint)=int((jul1-bdate)*86400.)
+        ireleaseend(numpoint)=int((jul2-bdate)*86400.)
+      else
+        ireleasestart(numpoint)=int((julm-bdate)*86400.)
+        ireleaseend(numpoint)=int((julm-bdate)*86400.)
+      endif
     else if (ldirect.eq.-1) then
       if ((jul1.lt.edate).or.(jul2.gt.bdate)) then
@@ -401 +410 @@
         stop
       endif
-      ireleasestart(numpoint)=int((jul1-bdate)*86400.)
-      ireleaseend(numpoint)=int((jul2-bdate)*86400.)
+      if (npart(numpoint).gt.num_min_discrete) then
+        ireleasestart(numpoint)=int((jul1-bdate)*86400.)
+        ireleaseend(numpoint)=int((jul2-bdate)*86400.)
+      else
+        ireleasestart(numpoint)=int((julm-bdate)*86400.)
+        ireleaseend(numpoint)=int((julm-bdate)*86400.)
+      endif
     endif
   endif