Index: trunk/src/readoutgrid.f90
===================================================================
--- trunk/src/readoutgrid.f90	(revision 20)
+++ trunk/src/readoutgrid.f90	(revision 27)
@@ -29,4 +29,6 @@
   !                                                                            *
   !     4 June 1996                                                            *
+  !     HSO, 1 July 2014
+  !     Added optional namelist input
   !                                                                            *
   !*****************************************************************************
@@ -53,32 +55,58 @@
   real,parameter :: eps=1.e-4
 
-
+  ! namelist variables
+  integer, parameter :: maxoutlev=500
+  integer :: readerror
+  real,allocatable, dimension (:) :: outheights
+
+  ! declare namelist
+  namelist /outgrid/ &
+    outlon0,outlat0, &
+    numxgrid,numygrid, &
+    dxout,dyout, &
+    outheights
+
+  ! allocate large array for reading input
+  allocate(outheights(maxoutlev),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate outheights'
+
+  ! helps identifying failed namelist input
+  dxout=-1.0
+  outheights=-1.0
 
   ! Open the OUTGRID file and read output grid specifications
   !**********************************************************
 
-  open(unitoutgrid,file=path(1)(1:length(1))//'OUTGRID',status='old', &
-       err=999)
-
-
-  call skplin(5,unitoutgrid)
-
-
-  ! 1.  Read horizontal grid specifications
-  !****************************************
-
-  call skplin(3,unitoutgrid)
-  read(unitoutgrid,'(4x,f11.4)') outlon0
-  call skplin(3,unitoutgrid)
-  read(unitoutgrid,'(4x,f11.4)') outlat0
-  call skplin(3,unitoutgrid)
-  read(unitoutgrid,'(4x,i5)') numxgrid
-  call skplin(3,unitoutgrid)
-  read(unitoutgrid,'(4x,i5)') numygrid
-  call skplin(3,unitoutgrid)
-  read(unitoutgrid,'(4x,f12.5)') dxout
-  call skplin(3,unitoutgrid)
-  read(unitoutgrid,'(4x,f12.5)') dyout
-
+  open(unitoutgrid,file=path(1)(1:length(1))//'OUTGRID',status='old',form='formatted',err=999)
+
+  ! try namelist input
+  read(unitoutgrid,outgrid,iostat=readerror)
+  close(unitoutgrid)
+
+  if ((dxout.le.0).or.(readerror.ne.0)) then
+
+    readerror=1
+
+    open(unitoutgrid,file=path(1)(1:length(1))//'OUTGRID',status='old',err=999)
+
+    call skplin(5,unitoutgrid)
+
+    ! 1.  Read horizontal grid specifications
+    !****************************************
+
+    call skplin(3,unitoutgrid)
+    read(unitoutgrid,'(4x,f11.4)') outlon0
+    call skplin(3,unitoutgrid)
+    read(unitoutgrid,'(4x,f11.4)') outlat0
+    call skplin(3,unitoutgrid)
+    read(unitoutgrid,'(4x,i5)') numxgrid
+    call skplin(3,unitoutgrid)
+    read(unitoutgrid,'(4x,i5)') numygrid
+    call skplin(3,unitoutgrid)
+    read(unitoutgrid,'(4x,f12.5)') dxout
+    call skplin(3,unitoutgrid)
+    read(unitoutgrid,'(4x,f12.5)') dyout
+
+  endif
 
   ! Check validity of output grid (shall be within model domain)
@@ -91,7 +119,5 @@
   if ((outlon0+eps.lt.xlon0).or.(outlat0+eps.lt.ylat0) &
        .or.(xr.gt.xr1+eps).or.(yr.gt.yr1+eps)) then
-    write(*,*) 'outlon0,outlat0:'
     write(*,*) outlon0,outlat0
-    write(*,*) 'xr1,yr1,xlon0,ylat0,xr,yr,dxout,dyout:'
     write(*,*) xr1,yr1,xlon0,ylat0,xr,yr,dxout,dyout
     write(*,*) ' #### FLEXPART MODEL ERROR! PART OF OUTPUT    ####'
@@ -104,6 +130,8 @@
   ! 2. Count Vertical levels of output grid
   !****************************************
-  j=0
-100   j=j+1
+
+  if (readerror.ne.0) then
+    j=0
+100 j=j+1
     do i=1,3
       read(unitoutgrid,*,end=99)
@@ -112,31 +140,58 @@
     if (outhelp.eq.0.) goto 99
     goto 100
-99   numzgrid=j-1
-
-    allocate(outheight(numzgrid) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
-    allocate(outheighthalf(numzgrid) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
-
-
-  rewind(unitoutgrid)
-  call skplin(29,unitoutgrid)
+99  numzgrid=j-1
+  else
+    do i=1,maxoutlev
+      if (outheights(i).lt.0) exit
+    end do
+    numzgrid=i-1
+  end if
+
+  allocate(outheight(numzgrid),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate outheight'
+  allocate(outheighthalf(numzgrid),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate outheighthalf'
 
   ! 2. Vertical levels of output grid
   !**********************************
 
-  j=0
-1000   j=j+1
-    do i=1,3
-      read(unitoutgrid,*,end=990)
-    end do
-    read(unitoutgrid,'(4x,f7.1)',end=990) outhelp
-    if (outhelp.eq.0.) goto 99
-    outheight(j)=outhelp
-    goto 1000
-990   numzgrid=j-1
-
+  if (readerror.ne.0) then
+
+    rewind(unitoutgrid)
+    call skplin(29,unitoutgrid)
+
+    do j=1,numzgrid
+      do i=1,3
+        read(unitoutgrid,*)
+      end do
+      read(unitoutgrid,'(4x,f7.1)') outhelp
+      outheight(j)=outhelp
+      outheights(j)=outhelp
+    end do
+    close(unitoutgrid)
+
+  else
+
+    do j=1,numzgrid
+      outheight(j)=outheights(j)
+    end do
+
+  endif
+
+  ! write outgrid file in namelist format to output directory if requested
+  if (nmlout.eqv..true.) then
+    ! reallocate outheights with actually required dimension for namelist writing
+    deallocate(outheights)
+    allocate(outheights(numzgrid),stat=stat)
+    if (stat.ne.0) write(*,*)'ERROR: could not allocate outheights'
+
+    do j=1,numzgrid
+      outheights(j)=outheight(j)
+    end do
+
+    open(unitoutgrid,file=path(2)(1:length(2))//'OUTGRID.namelist',err=1000)
+    write(unitoutgrid,nml=outgrid)
+    close(unitoutgrid)
+  endif
 
   ! Check whether vertical levels are specified in ascending order
@@ -160,31 +215,29 @@
   end do
 
-
   xoutshift=xlon0-outlon0
   youtshift=ylat0-outlat0
-  close(unitoutgrid)
-
-    allocate(oroout(0:numxgrid-1,0:numygrid-1) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
-    allocate(area(0:numxgrid-1,0:numygrid-1) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
-    allocate(volume(0:numxgrid-1,0:numygrid-1,numzgrid) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
-    allocate(areaeast(0:numxgrid-1,0:numygrid-1,numzgrid) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
-    allocate(areanorth(0:numxgrid-1,0:numygrid-1,numzgrid) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
+
+  allocate(oroout(0:numxgrid-1,0:numygrid-1),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate oroout'
+  allocate(area(0:numxgrid-1,0:numygrid-1),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate area'
+  allocate(volume(0:numxgrid-1,0:numygrid-1,numzgrid),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate volume'
+  allocate(areaeast(0:numxgrid-1,0:numygrid-1,numzgrid),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate areaeast'
+  allocate(areanorth(0:numxgrid-1,0:numygrid-1,numzgrid),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate areanorth'
   return
 
 
-999   write(*,*) ' #### FLEXPART MODEL ERROR! FILE "OUTGRID"    #### '
+999 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "OUTGRID"    #### '
   write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
-  write(*,*) ' #### xxx/flexpart/options                    #### '
+  write(*,'(a)') path(1)(1:length(1))
   stop
 
+1000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "OUTGRID"    #### '
+  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
+  write(*,'(a)') path(2)(1:length(2))
+  stop
+
 end subroutine readoutgrid
