Changeset 37 for branches/petra/src/readspecies.f90

Timestamp:

Apr 22, 2015, 3:42:35 PM (9 years ago)

Author:

pesei

Message:

Wet dep quick fix and other small changes. Wet depo quick fix not final yet.

File:

• branches/petra/src/readspecies.f90 (modified) (7 diffs)

branches/petra/src/readspecies.f90

@@ -20 +20 @@
 !**********************************************************************
 
-subroutine readspecies(id_spec,pos_spec)
+subroutine readspecies(id_spec,id_pos)
 
   !*****************************************************************************
   !                                                                            *
   !     This routine reads names and physical constants of chemical species/   *
-  !     radionuclides given in the parameter pos_spec                          *
+  !     radionuclides given in the parameter id_pos                          *
   !                                                                            *
   !   Author: A. Stohl                                                         *
@@ -32 +32 @@
   !
   !   Changes:                                                                 *
+  !                                                                            *
   !   N. Kristiansen, 31.01.2013: Including parameters for in-cloud scavenging *
   !                                                                            *
   !   HSO, 13 August 2013: added optional namelist input
   !   PS, 2/2015: access= -> position=
+  !   PS, 4/2015: quick fix for bug associated with end=22
   !                                                                            *
   !*****************************************************************************
@@ -60 +62 @@
   implicit none
 
-  integer :: i, pos_spec,j
+  integer :: i, id_pos,j
   integer :: idow,ihour,id_spec
   character(len=3) :: aspecnumb
@@ -69 +71 @@
   real :: pdsigma, pdryvel, pweightmolar, pohreact, pspec_ass, pkao
   real :: pweta_in, pwetb_in, pwetc_in, pwetd_in
-  integer :: readerror
+  integer :: ios
 
   ! declare namelist
@@ -100 +102 @@
   ! Open the SPECIES file and read species names and properties
   !************************************************************
-  specnum(pos_spec)=id_spec
-  write(aspecnumb,'(i3.3)') specnum(pos_spec)
+  specnum(id_pos)=id_spec
+  write(aspecnumb,'(i3.3)') specnum(id_pos)
   open(unitspecies,file=path(1)(1:length(1))//'SPECIES/SPECIES_'//aspecnumb,status='old',form='formatted',err=998)
-  !write(*,*) 'reading SPECIES',specnum(pos_spec)
+  !write(*,*) 'reading SPECIES',specnum(id_pos)
 
   ASSSPEC=.FALSE.
 
   ! try namelist input
-  read(unitspecies,species_params,iostat=readerror)
+  read(unitspecies,species_params,iostat=ios)
   close(unitspecies)
    
-  if ((pweightmolar.eq.-789.0).or.(readerror.ne.0)) then ! no namelist found
-
-    readerror=1
-
-    open(unitspecies,file=path(1)(1:length(1))//'SPECIES/SPECIES_'//aspecnumb,status='old',err=998)
+  if ((pweightmolar.eq.-789.0).or.(ios.ne.0)) then ! no namelist found
+
+    open(unitspecies,file=path(1) (1:length(1))//'SPECIES/SPECIES_'//aspecnumb, status='old', err=998)
 
     do i=1,6
       read(unitspecies,*)
     end do
-
-    read(unitspecies,'(a10)',end=22) species(pos_spec)
-  !  write(*,*) species(pos_spec)
-    read(unitspecies,'(f18.1)',end=22) decay(pos_spec)
-  !  write(*,*) decay(pos_spec)
-    read(unitspecies,'(e18.1)',end=22) weta(pos_spec)
-  !  write(*,*) weta(pos_spec)
-    read(unitspecies,'(f18.2)',end=22) wetb(pos_spec)
-  !  write(*,*) wetb(pos_spec)
+    read(unitspecies,'(a10)',end=20) species(id_pos)
+  !  write(*,*) species(id_pos)
+    read(unitspecies,'(f18.1)',end=20) decay(id_pos)
+  !  write(*,*) decay(id_pos)
+    read(unitspecies,'(e18.1)',end=20) weta(id_pos)
+  !  write(*,*) weta(id_pos)
+    read(unitspecies,'(f18.2)',end=20) wetb(id_pos)
+  !  write(*,*) wetb(id_pos)
 
   !*** NIK 31.01.2013: including in-cloud scavening parameters
-   read(unitspecies,'(e18.1)',end=22) weta_in(pos_spec)
-  !  write(*,*) weta_in(pos_spec)
-   read(unitspecies,'(f18.2)',end=22) wetb_in(pos_spec)
-  !  write(*,*) wetb_in(pos_spec)
-   read(unitspecies,'(f18.2)',end=22) wetc_in(pos_spec)
-  !  write(*,*) wetc_in(pos_spec)
-   read(unitspecies,'(f18.2)',end=22) wetd_in(pos_spec)
-  !  write(*,*) wetd_in(pos_spec)
-
-    read(unitspecies,'(f18.1)',end=22) reldiff(pos_spec)
-  !  write(*,*) reldiff(pos_spec)
-    read(unitspecies,'(e18.1)',end=22) henry(pos_spec)
-  !  write(*,*) henry(pos_spec)
-    read(unitspecies,'(f18.1)',end=22) f0(pos_spec)
-  !  write(*,*) f0(pos_spec)
-    read(unitspecies,'(e18.1)',end=22) density(pos_spec)
-  !  write(*,*) density(pos_spec)
-    read(unitspecies,'(e18.1)',end=22) dquer(pos_spec)
-  !  write(*,*) dquer(pos_spec)
-    read(unitspecies,'(e18.1)',end=22) dsigma(pos_spec)
-  !  write(*,*) dsigma(pos_spec)
-    read(unitspecies,'(f18.2)',end=22) dryvel(pos_spec)
-  !  write(*,*) dryvel(pos_spec)
-    read(unitspecies,'(f18.2)',end=22) weightmolar(pos_spec)
-  !  write(*,*) weightmolar(pos_spec)
-    read(unitspecies,'(e18.1)',end=22) ohreact(pos_spec)
-  !  write(*,*) ohreact(pos_spec)
-    read(unitspecies,'(i18)',end=22) spec_ass(pos_spec)
-  !  write(*,*) spec_ass(pos_spec)
-    read(unitspecies,'(f18.2)',end=22) kao(pos_spec)
-  !       write(*,*) kao(pos_spec)
-
-    pspecies=species(pos_spec)
-    pdecay=decay(pos_spec)
-    pweta=weta(pos_spec)
-    pwetb=wetb(pos_spec)
-    pweta_in=weta_in(pos_spec)
-    pwetb_in=wetb_in(pos_spec)
-    pwetc_in=wetc_in(pos_spec)
-    pwetd_in=wetd_in(pos_spec)
-    preldiff=reldiff(pos_spec)
-    phenry=henry(pos_spec)
-    pf0=f0(pos_spec)
-    pdensity=density(pos_spec)
-    pdquer=dquer(pos_spec)
-    pdsigma=dsigma(pos_spec)
-    pdryvel=dryvel(pos_spec)
-    pweightmolar=weightmolar(pos_spec)
-    pohreact=ohreact(pos_spec)
-    pspec_ass=spec_ass(pos_spec)
-    pkao=kao(pos_spec)
+   read(unitspecies,'(e18.1)',end=20) weta_in(id_pos)
+  !  write(*,*) weta_in(id_pos)
+   read(unitspecies,'(f18.2)',end=20) wetb_in(id_pos)
+  !  write(*,*) wetb_in(id_pos)
+   read(unitspecies,'(f18.2)',end=20) wetc_in(id_pos)
+  !  write(*,*) wetc_in(id_pos)
+   read(unitspecies,'(f18.2)',end=20) wetd_in(id_pos)
+  !  write(*,*) wetd_in(id_pos)
+
+    read(unitspecies,'(f18.1)',end=20) reldiff(id_pos)
+  !  write(*,*) reldiff(id_pos)
+    read(unitspecies,'(e18.1)',end=20) henry(id_pos)
+  !  write(*,*) henry(id_pos)
+    read(unitspecies,'(f18.1)',end=20) f0(id_pos)
+  !  write(*,*) f0(id_pos)
+    read(unitspecies,'(e18.1)',end=20) density(id_pos)
+  !  write(*,*) density(id_pos)
+    read(unitspecies,'(e18.1)',end=20) dquer(id_pos)
+  !  write(*,*) dquer(id_pos)
+    read(unitspecies,'(e18.1)',end=20) dsigma(id_pos)
+  !  write(*,*) dsigma(id_pos)
+    read(unitspecies,'(f18.2)',end=20) dryvel(id_pos)
+  !  write(*,*) dryvel(id_pos)
+    read(unitspecies,'(f18.2)',end=20) weightmolar(id_pos)
+  !  write(*,*) weightmolar(id_pos)
+    read(unitspecies,'(e18.1)',end=20) ohreact(id_pos)
+  !  write(*,*) ohreact(id_pos)
+    read(unitspecies,'(i18)',end=20) spec_ass(id_pos)
+  !  write(*,*) spec_ass(id_pos)
+    read(unitspecies,'(f18.2)',end=20) kao(id_pos)
+  !       write(*,*) kao(id_pos)
+
+    pspecies=species(id_pos)
+    pdecay=decay(id_pos)
+    pweta=weta(id_pos)
+    pwetb=wetb(id_pos)
+    pweta_in=weta_in(id_pos)
+    pwetb_in=wetb_in(id_pos)
+    pwetc_in=wetc_in(id_pos)
+    pwetd_in=wetd_in(id_pos)
+    preldiff=reldiff(id_pos)
+    phenry=henry(id_pos)
+    pf0=f0(id_pos)
+    pdensity=density(id_pos)
+    pdquer=dquer(id_pos)
+    pdsigma=dsigma(id_pos)
+    pdryvel=dryvel(id_pos)
+    pweightmolar=weightmolar(id_pos)
+    pohreact=ohreact(id_pos)
+    pspec_ass=spec_ass(id_pos)
+    pkao=kao(id_pos)
 
   else
 
-    species(pos_spec)=pspecies
-    decay(pos_spec)=pdecay
-    weta(pos_spec)=pweta
-    wetb(pos_spec)=pwetb
-    weta_in(pos_spec)=pweta_in
-    wetb_in(pos_spec)=pwetb_in
-    wetc_in(pos_spec)=pwetc_in
-    wetd_in(pos_spec)=pwetd_in
-    reldiff(pos_spec)=preldiff
-    henry(pos_spec)=phenry
-    f0(pos_spec)=pf0
-    density(pos_spec)=pdensity
-    dquer(pos_spec)=pdquer
-    dsigma(pos_spec)=pdsigma
-    dryvel(pos_spec)=pdryvel
-    weightmolar(pos_spec)=pweightmolar
-    ohreact(pos_spec)=pohreact
-    spec_ass(pos_spec)=pspec_ass
-    kao(pos_spec)=pkao
-
-  endif
-
-  i=pos_spec
-
-  if ((weta(pos_spec).gt.0).and.(henry(pos_spec).le.0)) then
-   if (dquer(pos_spec).le.0) goto 996 ! no particle, no henry set
-  endif
-
-  if (spec_ass(pos_spec).gt.0) then
+    species(id_pos)=pspecies
+    decay(id_pos)=pdecay
+    weta(id_pos)=pweta
+    wetb(id_pos)=pwetb
+    weta_in(id_pos)=pweta_in
+    wetb_in(id_pos)=pwetb_in
+    wetc_in(id_pos)=pwetc_in
+    wetd_in(id_pos)=pwetd_in
+    reldiff(id_pos)=preldiff
+    henry(id_pos)=phenry
+    f0(id_pos)=pf0
+    density(id_pos)=pdensity
+    dquer(id_pos)=pdquer
+    dsigma(id_pos)=pdsigma
+    dryvel(id_pos)=pdryvel
+    weightmolar(id_pos)=pweightmolar
+    ohreact(id_pos)=pohreact
+    spec_ass(id_pos)=pspec_ass
+    kao(id_pos)=pkao
+
+  endif
+
+  if ((weta(id_pos).gt.0).and.(henry(id_pos).le.0)) then
+   if (dquer(id_pos).le.0) goto 996 ! no particle, no henry set
+  endif
+
+  if (spec_ass(id_pos).gt.0) then
     spec_found=.FALSE.
-    do j=1,pos_spec-1
-      if (spec_ass(pos_spec).eq.specnum(j)) then
-        spec_ass(pos_spec)=j
+    do j=1,id_pos-1
+      if (spec_ass(id_pos).eq.specnum(j)) then
+        spec_ass(id_pos)=j
         spec_found=.TRUE.
         ASSSPEC=.TRUE.
@@ -227 +224 @@
   endif
 
-  if (dsigma(i).eq.1.) dsigma(i)=1.0001   ! avoid floating exception
-  if (dsigma(i).eq.0.) dsigma(i)=1.0001   ! avoid floating exception
-
-  if ((reldiff(i).gt.0.).and.(density(i).gt.0.)) then
+  if (dsigma(id_pos).eq.1.) dsigma(id_pos)=1.0001   ! avoid floating exception
+  if (dsigma(id_pos).eq.0.) dsigma(id_pos)=1.0001   ! avoid floating exception
+
+  if ((reldiff(id_pos).gt.0.).and.(density(id_pos).gt.0.)) then
     write(*,*) '#### FLEXPART MODEL ERROR! FILE "SPECIES"    ####'
     write(*,*) '#### IS CORRUPT. SPECIES CANNOT BE BOTH      ####'
@@ -246 +243 @@
 
   do j=1,24           ! initialize everything to no variation
-    area_hour(i,j)=1.
-    point_hour(i,j)=1.
+    area_hour(id_pos,j)=1.
+    point_hour(id_pos,j)=1.
   end do
   do j=1,7
-    area_dow(i,j)=1.
-    point_dow(i,j)=1.
+    area_dow(id_pos,j)=1.
+    point_dow(id_pos,j)=1.
   end do
 
-  if (readerror.ne.0) then ! text format input
-
-    read(unitspecies,*,end=22)
+  if (ios.ne.0) then ! text format input
+
+    read(unitspecies,*,iostat=ios)
+    if (ios .ne. 0) goto 22 ! ifort does not like err=
     do j=1,24     ! 24 hours, starting with 0-1 local time
-      read(unitspecies,*) ihour,area_hour(i,j),point_hour(i,j)
+      read(unitspecies,*) ihour,area_hour(id_pos,j),point_hour(id_pos,j)
     end do
     read(unitspecies,*)
     do j=1,7      ! 7 days of the week, starting with Monday
-      read(unitspecies,*) idow,area_dow(i,j),point_dow(i,j)
+      read(unitspecies,*) idow,area_dow(id_pos,j),point_dow(id_pos,j)
     end do