readoutgrid.f90 8.45 KB
Newer Older
Matthias Langer's avatar
 
Matthias Langer committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
!**********************************************************************
! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
!                                                                     *
! This file is part of FLEXPART.                                      *
!                                                                     *
! FLEXPART is free software: you can redistribute it and/or modify    *
! it under the terms of the GNU General Public License as published by*
! the Free Software Foundation, either version 3 of the License, or   *
! (at your option) any later version.                                 *
!                                                                     *
! FLEXPART is distributed in the hope that it will be useful,         *
! but WITHOUT ANY WARRANTY; without even the implied warranty of      *
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
! GNU General Public License for more details.                        *
!                                                                     *
! You should have received a copy of the GNU General Public License   *
! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
!**********************************************************************

subroutine readoutgrid

  !*****************************************************************************
  !                                                                            *
  !     This routine reads the user specifications for the output grid.        *
  !                                                                            *
  !     Author: A. Stohl                                                       *
  !                                                                            *
  !     4 June 1996                                                            *
pesei's avatar
pesei committed
31
32
33
  !     HSO, 1 July 2014: Add optional namelist input                          *
  !     PS, 6/2015-9/2018: read regular input with free format                 *
  !       simplify code and rename some variables                              *
Matthias Langer's avatar
 
Matthias Langer committed
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
  !                                                                            *
  !*****************************************************************************
  !                                                                            *
  ! Variables:                                                                 *
  ! dxout,dyout          grid distance                                         *
  ! numxgrid,numygrid,numzgrid    grid dimensions                              *
  ! outlon0,outlat0      lower left corner of grid                             *
  ! outheight(maxzgrid)  height levels of output grid [m]                      *
  !                                                                            *
  ! Constants:                                                                 *
  ! unitoutgrid          unit connected to file OUTGRID                        *
  !                                                                            *
  !*****************************************************************************

  use outg_mod
  use par_mod
  use com_mod

  implicit none

pesei's avatar
pesei committed
54
55
  integer :: i,kz,istat
  real :: xr,xr1,yr,yr1
Matthias Langer's avatar
 
Matthias Langer committed
56
57
  real,parameter :: eps=1.e-4

58
59
  ! namelist variables
  integer, parameter :: maxoutlev=500
pesei's avatar
pesei committed
60
61
62
  integer :: ios
  real,allocatable, dimension (:) :: outheights,outaux
  logical :: lnml
Matthias Langer's avatar
 
Matthias Langer committed
63

64
  ! declare namelist
pesei's avatar
pesei committed
65
  namelist /nml_outgrid/ &
66
67
68
    outlon0,outlat0, &
    numxgrid,numygrid, &
    dxout,dyout, &
pesei's avatar
pesei committed
69
    outheight
70

pesei's avatar
pesei committed
71
72
73
! allocate outheights for nml read with max dimension
  allocate(outheights(maxoutlev),outaux(maxoutlev),stat=istat)
  if (istat .ne. 0) write(*,*)'ERROR: could not allocate outheights'
Matthias Langer's avatar
 
Matthias Langer committed
74
75
76
77

  ! Open the OUTGRID file and read output grid specifications
  !**********************************************************

pesei's avatar
pesei committed
78
79
80
  outheight(:) = -999. ! initialise for later finding #valid levels
  open(unitoutgrid,file=trim(path(1))//'OUTGRID',status='old',&
    form='formatted',err=999)
Matthias Langer's avatar
 
Matthias Langer committed
81

82
  ! try namelist input
pesei's avatar
pesei committed
83
  read(unitoutgrid,nml_outgrid,iostat=ios)
84
  close(unitoutgrid)
Matthias Langer's avatar
 
Matthias Langer committed
85

pesei's avatar
pesei committed
86
87
88
  if (ios .eq. 0) then ! namelist works

    lnml = .true.
Matthias Langer's avatar
 
Matthias Langer committed
89

pesei's avatar
pesei committed
90
  else ! read as regular text
Matthias Langer's avatar
 
Matthias Langer committed
91

pesei's avatar
pesei committed
92
    lnml = .false.
Matthias Langer's avatar
 
Matthias Langer committed
93

pesei's avatar
pesei committed
94
    open(unitoutgrid,file=trim(path(1))//'OUTGRID',status='old',err=999)
95
    call skplin(5,unitoutgrid)
Matthias Langer's avatar
 
Matthias Langer committed
96

pesei's avatar
pesei committed
97
   ! Read horizontal grid specifications
98
99
100
    !****************************************

    call skplin(3,unitoutgrid)
pesei's avatar
pesei committed
101
    read(unitoutgrid,*) outlon0
102
    call skplin(3,unitoutgrid)
pesei's avatar
pesei committed
103
    read(unitoutgrid,*) outlat0
104
    call skplin(3,unitoutgrid)
pesei's avatar
pesei committed
105
    read(unitoutgrid,*) numxgrid
106
    call skplin(3,unitoutgrid)
pesei's avatar
pesei committed
107
    read(unitoutgrid,*) numygrid
108
    call skplin(3,unitoutgrid)
pesei's avatar
pesei committed
109
    read(unitoutgrid,*) dxout
110
    call skplin(3,unitoutgrid)
pesei's avatar
pesei committed
111
    read(unitoutgrid,*) dyout
112

pesei's avatar
pesei committed
113
  endif ! read OUTGRID file
Matthias Langer's avatar
 
Matthias Langer committed
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128

  ! Check validity of output grid (shall be within model domain)
  !*************************************************************

  xr=outlon0+real(numxgrid)*dxout
  yr=outlat0+real(numygrid)*dyout
  xr1=xlon0+real(nxmin1)*dx
  yr1=ylat0+real(nymin1)*dy
  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(*,*) xr1,yr1,xlon0,ylat0,xr,yr,dxout,dyout
    write(*,*) ' #### FLEXPART MODEL ERROR! PART OF OUTPUT    ####'
    write(*,*) ' #### GRID IS OUTSIDE MODEL DOMAIN. CHANGE    ####'
    write(*,*) ' #### FILE OUTGRID IN DIRECTORY               ####'
pesei's avatar
pesei committed
129
    write(*,'(a)') trim(path(1))
Matthias Langer's avatar
 
Matthias Langer committed
130
131
132
    stop
  endif

pesei's avatar
pesei committed
133
134
! Read (if .not. lmnl) and count vertical levels of output grid 
!**************************************************************
Matthias Langer's avatar
 
Matthias Langer committed
135

pesei's avatar
pesei committed
136
137
138
139
140
141
142
143
144
  do kz = 1,maxoutlev
    if (lnml) then ! we have read them already
      if (outheight(kz) .lt. 0.) exit ! 1st nondefined level
    else
      call skplin(3,unitoutgrid)
      read(unitoutgrid,*,end=10) outheight(kz)
    endif
  end do
10 continue  
145

pesei's avatar
pesei committed
146
  numzgrid = kz - 1 ! number of outgrid levels
147

pesei's avatar
pesei committed
148
149
! allocate the required length only, shuffle data  
  outaux = outheight ! shuffle
150

pesei's avatar
pesei committed
151
152
153
154
155
156
  deallocate(outheights)
  allocate(outheight(numzgrid),outheighthalf(numzgrid),stat=istat)
  if (istat .ne. 0) then
    write(*,*) 'ERROR: could not allocate outheight and outheighthalf'
    stop 'readoutgrid error'
  endif
157

pesei's avatar
pesei committed
158
159
  outheight=outaux(1:numzgrid) ! shuffle back
  deallocate (outaux) 
160

pesei's avatar
pesei committed
161
162
163
164
! write outgrid file in namelist format to output directory if requested
  if (nmlout) then
    open(unitoutgrid,file=trim(path(2))//'OUTGRID.namelist',err=1000)
    write(unitoutgrid,nml=nml_outgrid)
165
166
    close(unitoutgrid)
  endif
Matthias Langer's avatar
 
Matthias Langer committed
167
168
169
170

  ! Check whether vertical levels are specified in ascending order
  !***************************************************************

pesei's avatar
pesei committed
171
172
  do kz=2,numzgrid
    if (outheight(kz) .le. outheight(kz-1)) goto 998
Matthias Langer's avatar
 
Matthias Langer committed
173
174
175
176
177
  end do

  ! Determine the half levels, i.e. middle levels of the output grid
  !*****************************************************************

pesei's avatar
pesei committed
178
179
180
  outheighthalf(1) = 0.5*outheight(1)
  do kz = 2,numzgrid
    outheighthalf(kz) = 0.5*(outheight(kz-1)+outheight(kz))
Matthias Langer's avatar
 
Matthias Langer committed
181
182
183
184
185
  end do

  xoutshift=xlon0-outlon0
  youtshift=ylat0-outlat0

pesei's avatar
pesei committed
186
187
188
189
190
191
192
193
194
195
196
  allocate(oroout(0:numxgrid-1,0:numygrid-1),stat=istat)
  if (istat .ne. 0) write(*,*)'ERROR: could not allocate oroout'
  allocate(area(0:numxgrid-1,0:numygrid-1),stat=istat)
  if (istat .ne. 0) write(*,*)'ERROR: could not allocate area'
  allocate(volume(0:numxgrid-1,0:numygrid-1,numzgrid),stat=istat)
  if (istat .ne. 0) write(*,*)'ERROR: could not allocate volume'
  allocate(areaeast(0:numxgrid-1,0:numygrid-1,numzgrid),stat=istat)
  if (istat .ne. 0) write(*,*)'ERROR: could not allocate areaeast'
  allocate(areanorth(0:numxgrid-1,0:numygrid-1,numzgrid),stat=istat)
  if (istat .ne. 0) write(*,*)'ERROR: could not allocate areanorth'
  
Matthias Langer's avatar
 
Matthias Langer committed
197
198
  return

pesei's avatar
pesei committed
199
200
201
202
203
204
998 continue
  write(*,*) ' #### FLEXPART MODEL ERROR! YOUR SPECIFICATION#### '
  write(*,*) ' #### OF OUTPUT LEVELS NOT INCREASING AT LEVEL#### '
  write(*,*) ' #### ',kz,'                                  #### '
  write(*,*) ' #### PLEASE MAKE CHANGES IN FILE OUTGRID.    #### '
  STOP 'readoutgrid error'
Matthias Langer's avatar
 
Matthias Langer committed
205

206
207
999 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "OUTGRID"    #### '
  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
pesei's avatar
pesei committed
208
  write(*,'(a)') trim(path(1))
209
210
211
  stop

1000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "OUTGRID"    #### '
Matthias Langer's avatar
 
Matthias Langer committed
212
  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
pesei's avatar
pesei committed
213
  write(*,'(a)') trim(path(2))
Matthias Langer's avatar
 
Matthias Langer committed
214
215
216
  stop

end subroutine readoutgrid