Changeset 10 for branches/flexpart91_hasod/src_parallel/readreleases.f90

Timestamp:

Aug 15, 2013, 3:23:48 PM (11 years ago)

Author:

hasod

Message:

ADD: namelist input implemented for all common input files

File:

• branches/flexpart91_hasod/src_parallel/readreleases.f90 (modified) (8 diffs)

branches/flexpart91_hasod/src_parallel/readreleases.f90

@@ -33 +33 @@
   !     Update: 29 January 2001                                                *
   !     Release altitude can be either in magl or masl                         *
+  !     HSO, 12 August 2013
+  !     Added optional namelist input
   !                                                                            *
   !*****************************************************************************
@@ -67 +69 @@
   implicit none
 
-  integer :: numpartmax,i,j,id1,it1,id2,it2,specnum_rel,idum,stat
+  integer :: numpartmax,i,j,id1,it1,id2,it2,idum,stat
+  integer :: specnum_rel(maxspec)
   real :: vsh(ni),fracth(ni),schmih(ni),releaserate,xdum,cun
   real(kind=dp) :: jul1,jul2,juldate
@@ -73 +76 @@
   logical :: old
 
+  ! help variables for namelist reading
+  integer :: numpoints, parts, readerror
+  integer*2 :: zkind
+  integer :: idate1, itime1, idate2, itime2
+  real :: lon1,lon2,lat1,lat2,z1,z2,mass(nspec)
+  character*40 :: comment
+
+  ! declare namelists
+  namelist /releases_ctrl/ &
+    nspec, &
+    specnum_rel
+
+  namelist /release/ &
+    idate1, itime1, &
+    idate2, itime2, &
+    lon1, lon2, &
+    lat1, lat2, &
+    z1, z2, &
+    zkind, &
+    mass, &
+    parts, &
+    comment
+
+  numpoint=0
+
+  ! presetting namelist releases_ctrl
+  nspec = -1  ! use negative value to determine failed namelist input
+  specnum_rel(1) = 0
+
   !sec, read release to find how many releasepoints should be allocated
-
   open(unitreleases,file=path(1)(1:length(1))//'RELEASES',status='old', &
        err=999)
+
+  ! check if namelist input provided
+  read(unitreleases,releases_ctrl,iostat=readerror)
+
+  if ((readerror.ne.0).or.(nspec.lt.0)) then
+
+  ! no namelist format, close file and allow reopening in old format
+  rewind(unitreleases)
+
+  readerror=1 ! indicates old format
 
   ! Check the format of the RELEASES file (either in free format,
@@ -93 +134 @@
   rewind(unitreleases)
 
-
   ! Skip first 11 lines (file header)
   !**********************************
 
   call skplin(11,unitreleases)
-
-
   read(unitreleases,*,err=998) nspec
   if (old) call skplin(2,unitreleases)
   do i=1,nspec
-    read(unitreleases,*,err=998) specnum_rel
+    read(unitreleases,*,err=998) specnum_rel(1)
     if (old) call skplin(2,unitreleases)
   end do
 
-  numpoint=0
 100   numpoint=numpoint+1
   read(unitreleases,*,end=25)
@@ -142 +179 @@
 25   numpoint=numpoint-1
 
+  else 
+    readerror=0
+    do while (readerror.eq.0) 
+      idate1=-1
+      read(unitreleases,release,iostat=readerror)
+      if ((idate1.lt.0).or.(readerror.ne.0)) then
+        readerror=1
+      else
+        numpoint=numpoint+1
+      endif
+    end do
+    readerror=0
+  endif ! if namelist input
+
+  rewind(unitreleases)
+
   !allocate memory for numpoint releaspoint
-    allocate(ireleasestart(numpoint) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
-    allocate(ireleaseend(numpoint) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
-    allocate(xpoint1(numpoint) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
-    allocate(xpoint2(numpoint) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
-    allocate(ypoint1(numpoint) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
-    allocate(ypoint2(numpoint) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
-    allocate(zpoint1(numpoint) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
-    allocate(zpoint2(numpoint) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
-    allocate(kindz(numpoint) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
-    allocate(xmass(numpoint,maxspec) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
-    allocate(rho_rel(numpoint) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
-    allocate(npart(numpoint) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
-    allocate(xmasssave(numpoint) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
-
-   write (*,*) ' Releasepoints allocated: ', numpoint
-
-   do i=1,numpoint
-     xmasssave(i)=0.
-   end do
+  allocate(ireleasestart(numpoint),stat=stat)
+       
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
+  allocate(ireleaseend(numpoint),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
+  allocate(xpoint1(numpoint),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
+  allocate(xpoint2(numpoint),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
+  allocate(ypoint1(numpoint),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
+  allocate(ypoint2(numpoint),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
+  allocate(zpoint1(numpoint),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
+  allocate(zpoint2(numpoint),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
+  allocate(kindz(numpoint),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
+  allocate(xmass(numpoint,maxspec),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
+  allocate(rho_rel(numpoint),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
+  allocate(npart(numpoint),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
+  allocate(xmasssave(numpoint),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
+
+  write (*,*) numpoint,' releasepoints allocated.'
+
+  do i=1,numpoint
+    xmasssave(i)=0.
+  end do
 
   !now save the information
@@ -198 +239 @@
   end do
 
-  rewind(unitreleases)
-
-
+  if (readerror.ne.0) then
   ! Skip first 11 lines (file header)
   !**********************************
 
-  call skplin(11,unitreleases)
-
+    call skplin(11,unitreleases)
 
   ! Assign species-specific parameters needed for physical processes
   !*************************************************************************
 
-  read(unitreleases,*,err=998) nspec
-  if (nspec.gt.maxspec) goto 994
-  if (old) call skplin(2,unitreleases)
+    read(unitreleases,*,err=998) nspec
+    if (nspec.gt.maxspec) goto 994
+  if (old) call skplin(2,unitreleases)
+  endif
   do i=1,nspec
-    read(unitreleases,*,err=998) specnum_rel
-    if (old) call skplin(2,unitreleases)
-    call readspecies(specnum_rel,i)
+    if (readerror.ne.0) then
+      read(unitreleases,*,err=998) specnum_rel(1)
+      if (old) call skplin(2,unitreleases)
+      call readspecies(specnum_rel(1),i)
+    else
+      call readspecies(specnum_rel(i),i)
+    endif
 
   ! For backward runs, only 1 species is allowed
@@ -301 +344 @@
   ! Read specifications for each release point
   !*******************************************
-
+  numpoints=numpoint
   numpoint=0
   numpartmax=0
   releaserate=0.
 1000   numpoint=numpoint+1
+
+  if (readerror.eq.0) then ! reading namelist format
+
+  if (numpoint.gt.numpoints) goto 250
+  zkind = 1
+  mass = 0
+  parts = 0
+  comment = ' '
+  read(unitreleases,release,iostat=readerror)
+  id1=idate1
+  it1=itime1
+  id2=idate2
+  it2=itime2
+  xpoint1(numpoint)=lon1
+  xpoint2(numpoint)=lon2
+  ypoint1(numpoint)=lat1
+  ypoint2(numpoint)=lat2
+  zpoint1(numpoint)=z1
+  zpoint2(numpoint)=z2
+  kindz(numpoint)=zkind
+  do i=1,nspec
+    xmass(numpoint,i)=mass(i)
+  end do
+  npart(numpoint)=parts
+  compoint(min(1001,numpoint))=comment
+
+  else
+
   read(unitreleases,*,end=250)
   read(unitreleases,*,err=998,end=250) id1,it1
@@ -347 +418 @@
   endif
 
+  endif ! if namelist format
 
   ! If a release point contains no particles, stop and issue error message