source: trunk/src/outgrid_init.f90 @ 24

Last change on this file since 24 was 20, checked in by igpis, 10 years ago

move version 9.1.8 form branches to trunk. Contributions from HSO, saeck, pesei, NIK, RT, XKF, IP and others
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 nested and global output (x,y): ', &
228       max(numxgrid,numxgridn),max(numygrid,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