source: trunk/src/outgrid_init.f90 @ 28

Last change on this file since 28 was 27, checked in by hasod, 10 years ago
  • Implemented optional namelist input for COMMAND, RELEASES, SPECIES, AGECLASSES,OUTGRID,OUTGRID_NEST,RECEPTORS
  • Implemented com_mod switch nmlout to write input files as namelist to the output directory (.true. by default)
  • Proposed updated startup and runtime output (may change back to previous info if desired)
File size: 11.9 KB
Line 
1!**********************************************************************
2! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
3! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
4! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
5!                                                                     *
6! This file is part of FLEXPART.                                      *
7!                                                                     *
8! FLEXPART is free software: you can redistribute it and/or modify    *
9! it under the terms of the GNU General Public License as published by*
10! the Free Software Foundation, either version 3 of the License, or   *
11! (at your option) any later version.                                 *
12!                                                                     *
13! FLEXPART is distributed in the hope that it will be useful,         *
14! but WITHOUT ANY WARRANTY; without even the implied warranty of      *
15! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
16! GNU General Public License for more details.                        *
17!                                                                     *
18! You should have received a copy of the GNU General Public License   *
19! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
20!**********************************************************************
21
22subroutine outgrid_init
23  !
24  !*****************************************************************************
25  !                                                                            *
26  !  This routine initializes the output grids                                 *
27  !                                                                            *
28  !     Author: A. Stohl                                                       *
29  !                                                                            *
30  !     7 August 2002                                                          *
31  !                                                                            *
32  !*****************************************************************************
33  !                                                                            *
34  ! Variables:                                                                 *
35  !                                                                            *
36  ! area               surface area of all output grid cells                   *
37  ! areaeast           eastward facing wall area of all output grid cells      *
38  ! areanorth          northward facing wall area of all output grid cells     *
39  ! volume             volumes of all output grid cells                        *
40  !                                                                            *
41  !*****************************************************************************
42
43  use flux_mod
44  use oh_mod
45  use unc_mod
46  use outg_mod
47  use par_mod
48  use com_mod
49
50  implicit none
51
52  integer :: ix,jy,kz,i,nage,l,iix,jjy,ixp,jyp,i1,j1,j,ngrid
53  integer :: ks,kp,stat
54  real :: ylat,gridarea,ylatp,ylatm,hzone,cosfactm,cosfactp
55  real :: xlon,xl,yl,ddx,ddy,rddx,rddy,p1,p2,p3,p4,xtn,ytn,oroh
56  real,parameter :: eps=nxmax/3.e5
57
58
59  ! Compute surface area and volume of each grid cell: area, volume;
60  ! and the areas of the northward and eastward facing walls: areaeast, areanorth
61  !***********************************************************************
62  do jy=0,numygrid-1
63    ylat=outlat0+(real(jy)+0.5)*dyout
64    ylatp=ylat+0.5*dyout
65    ylatm=ylat-0.5*dyout
66    if ((ylatm.lt.0).and.(ylatp.gt.0.)) then
67      hzone=dyout*r_earth*pi180
68    else
69
70  ! Calculate area of grid cell with formula M=2*pi*R*h*dx/360,
71  ! see Netz, Formeln der Mathematik, 5. Auflage (1983), p.90
72  !************************************************************
73
74      cosfactp=cos(ylatp*pi180)
75      cosfactm=cos(ylatm*pi180)
76      if (cosfactp.lt.cosfactm) then
77        hzone=sqrt(1-cosfactp**2)- &
78             sqrt(1-cosfactm**2)
79        hzone=hzone*r_earth
80      else
81        hzone=sqrt(1-cosfactm**2)- &
82             sqrt(1-cosfactp**2)
83        hzone=hzone*r_earth
84      endif
85    endif
86
87  ! Surface are of a grid cell at a latitude ylat
88  !**********************************************
89
90    gridarea=2.*pi*r_earth*hzone*dxout/360.
91
92    do ix=0,numxgrid-1
93      area(ix,jy)=gridarea
94
95  ! Volume = area x box height
96  !***************************
97
98      volume(ix,jy,1)=area(ix,jy)*outheight(1)
99      areaeast(ix,jy,1)=dyout*r_earth*pi180*outheight(1)
100      areanorth(ix,jy,1)=cos(ylat*pi180)*dxout*r_earth*pi180* &
101           outheight(1)
102      do kz=2,numzgrid
103        areaeast(ix,jy,kz)=dyout*r_earth*pi180* &
104             (outheight(kz)-outheight(kz-1))
105        areanorth(ix,jy,kz)=cos(ylat*pi180)*dxout*r_earth*pi180* &
106             (outheight(kz)-outheight(kz-1))
107        volume(ix,jy,kz)=area(ix,jy)*(outheight(kz)-outheight(kz-1))
108      end do
109    end do
110  end do
111
112
113
114
115  !******************************************************************
116  ! Determine average height of model topography in output grid cells
117  !******************************************************************
118
119  ! Loop over all output grid cells
120  !********************************
121
122  do jjy=0,numygrid-1
123    do iix=0,numxgrid-1
124      oroh=0.
125
126  ! Take 100 samples of the topography in every grid cell
127  !******************************************************
128
129      do j1=1,10
130        ylat=outlat0+(real(jjy)+real(j1)/10.-0.05)*dyout
131        yl=(ylat-ylat0)/dy
132        do i1=1,10
133          xlon=outlon0+(real(iix)+real(i1)/10.-0.05)*dxout
134          xl=(xlon-xlon0)/dx
135
136  ! Determine the nest we are in
137  !*****************************
138
139          ngrid=0
140          do j=numbnests,1,-1
141            if ((xl.gt.xln(j)+eps).and.(xl.lt.xrn(j)-eps).and. &
142                 (yl.gt.yln(j)+eps).and.(yl.lt.yrn(j)-eps)) then
143              ngrid=j
144              goto 43
145            endif
146          end do
14743        continue
148
149  ! Determine (nested) grid coordinates and auxiliary parameters used for interpolation
150  !*****************************************************************************
151
152          if (ngrid.gt.0) then
153            xtn=(xl-xln(ngrid))*xresoln(ngrid)
154            ytn=(yl-yln(ngrid))*yresoln(ngrid)
155            ix=int(xtn)
156            jy=int(ytn)
157            ddy=ytn-real(jy)
158            ddx=xtn-real(ix)
159          else
160            ix=int(xl)
161            jy=int(yl)
162            ddy=yl-real(jy)
163            ddx=xl-real(ix)
164          endif
165          ixp=ix+1
166          jyp=jy+1
167          rddx=1.-ddx
168          rddy=1.-ddy
169          p1=rddx*rddy
170          p2=ddx*rddy
171          p3=rddx*ddy
172          p4=ddx*ddy
173
174          if (ngrid.gt.0) then
175            oroh=oroh+p1*oron(ix ,jy ,ngrid) &
176                 + p2*oron(ixp,jy ,ngrid) &
177                 + p3*oron(ix ,jyp,ngrid) &
178                 + p4*oron(ixp,jyp,ngrid)
179          else
180            oroh=oroh+p1*oro(ix ,jy) &
181                 + p2*oro(ixp,jy) &
182                 + p3*oro(ix ,jyp) &
183                 + p4*oro(ixp,jyp)
184          endif
185        end do
186      end do
187
188  ! Divide by the number of samples taken
189  !**************************************
190
191      oroout(iix,jjy)=oroh/100.
192    end do
193  end do
194
195  ! if necessary allocate flux fields
196  if (iflux.eq.1) then
197    allocate(flux(6,0:numxgrid-1,0:numygrid-1,numzgrid, &
198         1:nspec,1:maxpointspec_act,1:nageclass),stat=stat)
199    if (stat.ne.0) write(*,*)'ERROR: could not allocate flux array '
200  endif
201
202  !write (*,*) 'allocating: in a sec',OHREA
203  if (OHREA.eqv..TRUE.) then
204  !   write (*,*) 'allocating: ',maxxOH,maxyOH,maxzOH
205    allocate(OH_field(12,0:maxxOH-1,0:maxyOH-1,maxzOH) &
206         ,stat=stat)
207    if (stat.ne.0) write(*,*)'ERROR: could not allocate OH array '
208    allocate(OH_field_height(7) &
209         ,stat=stat)
210    if (stat.ne.0) write(*,*)'ERROR: could not allocate OH array '
211  endif
212  ! gridunc,griduncn        uncertainty of outputted concentrations
213  allocate(gridunc(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec, &
214       maxpointspec_act,nclassunc,maxageclass),stat=stat)
215    if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
216  if (ldirect.gt.0) then
217  allocate(wetgridunc(0:numxgrid-1,0:numygrid-1,maxspec, &
218       maxpointspec_act,nclassunc,maxageclass),stat=stat)
219    if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
220  allocate(drygridunc(0:numxgrid-1,0:numygrid-1,maxspec, &
221       maxpointspec_act,nclassunc,maxageclass),stat=stat)
222    if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
223  endif
224  !write (*,*) 'Dimensions for fields', numxgrid,numygrid, &
225  !     maxspec,maxpointspec_act,nclassunc,maxageclass
226
227  write (*,*) 'Allocating fields for global output (x,y): ', numxgrid,numygrid
228  write (*,*) 'Allocating fields for nested output (x,y): ', numxgridn,numygridn
229
230  ! allocate fields for concoutput with maximum dimension of outgrid
231  ! and outgrid_nest
232
233  allocate(gridsigma(0:max(numxgrid,numxgridn)-1, &
234       0:max(numygrid,numygridn)-1,numzgrid),stat=stat)
235    if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
236  allocate(grid(0:max(numxgrid,numxgridn)-1, &
237       0:max(numygrid,numygridn)-1,numzgrid),stat=stat)
238    if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
239  allocate(densityoutgrid(0:max(numxgrid,numxgridn)-1, &
240       0:max(numygrid,numygridn)-1,numzgrid),stat=stat)
241    if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
242
243  allocate(factor3d(0:max(numxgrid,numxgridn)-1, &
244       0:max(numygrid,numygridn)-1,numzgrid),stat=stat)
245    if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
246  allocate(sparse_dump_r(max(numxgrid,numxgridn)* &
247       max(numygrid,numygridn)*numzgrid),stat=stat)
248    if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
249
250   allocate(sparse_dump_u(max(numxgrid,numxgridn)* &
251       max(numygrid,numygridn)*numzgrid),stat=stat)
252        if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
253
254  allocate(sparse_dump_i(max(numxgrid,numxgridn)* &
255       max(numygrid,numygridn)*numzgrid),stat=stat)
256    if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
257
258  ! deposition fields are only allocated for forward runs
259  if (ldirect.gt.0) then
260     allocate(wetgridsigma(0:max(numxgrid,numxgridn)-1, &
261          0:max(numygrid,numygridn)-1),stat=stat)
262     if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
263     allocate(drygridsigma(0:max(numxgrid,numxgridn)-1, &
264          0:max(numygrid,numygridn)-1),stat=stat)
265     if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
266     allocate(wetgrid(0:max(numxgrid,numxgridn)-1, &
267          0:max(numygrid,numygridn)-1),stat=stat)
268     if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
269     allocate(drygrid(0:max(numxgrid,numxgridn)-1, &
270          0:max(numygrid,numygridn)-1),stat=stat)
271     if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
272  endif
273
274  ! Initial condition field
275
276  if (linit_cond.gt.0) then
277    allocate(init_cond(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec, &
278         maxpointspec_act),stat=stat)
279    if (stat.ne.0) write(*,*)'ERROR: could not allocate init_cond'
280  endif
281
282  !************************
283  ! Initialize output grids
284  !************************
285
286  do ks=1,nspec
287  do kp=1,maxpointspec_act
288    do i=1,numreceptor
289  ! Receptor points
290      creceptor(i,ks)=0.
291    end do
292    do nage=1,nageclass
293      do jy=0,numygrid-1
294        do ix=0,numxgrid-1
295          do l=1,nclassunc
296  ! Deposition fields
297            if (ldirect.gt.0) then
298            wetgridunc(ix,jy,ks,kp,l,nage)=0.
299            drygridunc(ix,jy,ks,kp,l,nage)=0.
300            endif
301            do kz=1,numzgrid
302              if (iflux.eq.1) then
303  ! Flux fields
304                 do i=1,5
305                   flux(i,ix,jy,kz,ks,kp,nage)=0.
306                 end do
307              endif
308  ! Initial condition field
309              if ((l.eq.1).and.(nage.eq.1).and.(linit_cond.gt.0)) &
310                   init_cond(ix,jy,kz,ks,kp)=0.
311  ! Concentration fields
312              gridunc(ix,jy,kz,ks,kp,l,nage)=0.
313            end do
314          end do
315        end do
316      end do
317    end do
318  end do
319  end do
320
321
322
323end subroutine outgrid_init
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG