source: flexpart.git/src/outgrid_init.f90 @ 02095e3

devrelease-10univie
Last change on this file since 02095e3 was 6ecb30a, checked in by Espen Sollum ATMOS <eso@…>, 2 years ago

Merged changes from CTBTO project

  • Property mode set to 100644
File size: 13.1 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!                                                                     *
22! DJM - 2017-05-09 - added #ifdef USE_MPIINPLACE cpp directive to     *
23! enable allocation of a gridunc0 array if required by MPI code in    *
24! mpi_mod.f90                                                         *
25!                                                                     *
26!**********************************************************************
27
28
29subroutine outgrid_init
30  !
31  !*****************************************************************************
32  !                                                                            *
33  !  This routine initializes the output grids                                 *
34  !                                                                            *
35  !     Author: A. Stohl                                                       *
36  !                                                                            *
37  !     7 August 2002                                                          *
38  !                                                                            *
39  !*****************************************************************************
40  !                                                                            *
41  ! Variables:                                                                 *
42  !                                                                            *
43  ! area               surface area of all output grid cells                   *
44  ! areaeast           eastward facing wall area of all output grid cells      *
45  ! areanorth          northward facing wall area of all output grid cells     *
46  ! volume             volumes of all output grid cells                        *
47  !                                                                            *
48  !*****************************************************************************
49
50  use flux_mod
51  use oh_mod
52  use unc_mod
53  use outg_mod
54  use par_mod
55  use com_mod
56
57  implicit none
58
59  integer :: ix,jy,kz,i,nage,l,iix,jjy,ixp,jyp,i1,j1,j,ngrid
60  integer :: ks,kp,stat
61  real :: ylat,gridarea,ylatp,ylatm,hzone,cosfactm,cosfactp
62  real :: xlon,xl,yl,ddx,ddy,rddx,rddy,p1,p2,p3,p4,xtn,ytn,oroh
63  real,parameter :: eps=nxmax/3.e5
64
65
66  ! Compute surface area and volume of each grid cell: area, volume;
67  ! and the areas of the northward and eastward facing walls: areaeast, areanorth
68  !***********************************************************************
69  do jy=0,numygrid-1
70    ylat=outlat0+(real(jy)+0.5)*dyout
71    ylatp=ylat+0.5*dyout
72    ylatm=ylat-0.5*dyout
73    if ((ylatm.lt.0).and.(ylatp.gt.0.)) then
74      hzone=dyout*r_earth*pi180
75    else
76
77  ! Calculate area of grid cell with formula M=2*pi*R*h*dx/360,
78  ! see Netz, Formeln der Mathematik, 5. Auflage (1983), p.90
79  !************************************************************
80
81      cosfactp=cos(ylatp*pi180)
82      cosfactm=cos(ylatm*pi180)
83      if (cosfactp.lt.cosfactm) then
84        hzone=sqrt(1-cosfactp**2)- &
85             sqrt(1-cosfactm**2)
86        hzone=hzone*r_earth
87      else
88        hzone=sqrt(1-cosfactm**2)- &
89             sqrt(1-cosfactp**2)
90        hzone=hzone*r_earth
91      endif
92    endif
93
94  ! Surface are of a grid cell at a latitude ylat
95  !**********************************************
96
97    gridarea=2.*pi*r_earth*hzone*dxout/360.
98
99    do ix=0,numxgrid-1
100      area(ix,jy)=gridarea
101
102  ! Volume = area x box height
103  !***************************
104
105      volume(ix,jy,1)=area(ix,jy)*outheight(1)
106      areaeast(ix,jy,1)=dyout*r_earth*pi180*outheight(1)
107      areanorth(ix,jy,1)=cos(ylat*pi180)*dxout*r_earth*pi180* &
108           outheight(1)
109      do kz=2,numzgrid
110        areaeast(ix,jy,kz)=dyout*r_earth*pi180* &
111             (outheight(kz)-outheight(kz-1))
112        areanorth(ix,jy,kz)=cos(ylat*pi180)*dxout*r_earth*pi180* &
113             (outheight(kz)-outheight(kz-1))
114        volume(ix,jy,kz)=area(ix,jy)*(outheight(kz)-outheight(kz-1))
115      end do
116    end do
117  end do
118
119
120
121
122  !******************************************************************
123  ! Determine average height of model topography in output grid cells
124  !******************************************************************
125
126  ! Loop over all output grid cells
127  !********************************
128
129  do jjy=0,numygrid-1
130    do iix=0,numxgrid-1
131      oroh=0.
132
133  ! Take 100 samples of the topography in every grid cell
134  !******************************************************
135
136      do j1=1,10
137        ylat=outlat0+(real(jjy)+real(j1)/10.-0.05)*dyout
138        yl=(ylat-ylat0)/dy
139        do i1=1,10
140          xlon=outlon0+(real(iix)+real(i1)/10.-0.05)*dxout
141          xl=(xlon-xlon0)/dx
142
143  ! Determine the nest we are in
144  !*****************************
145
146          ngrid=0
147          do j=numbnests,1,-1
148            if ((xl.gt.xln(j)+eps).and.(xl.lt.xrn(j)-eps).and. &
149                 (yl.gt.yln(j)+eps).and.(yl.lt.yrn(j)-eps)) then
150              ngrid=j
151              goto 43
152            endif
153          end do
15443        continue
155
156  ! Determine (nested) grid coordinates and auxiliary parameters used for interpolation
157  !*****************************************************************************
158
159          if (ngrid.gt.0) then
160            xtn=(xl-xln(ngrid))*xresoln(ngrid)
161            ytn=(yl-yln(ngrid))*yresoln(ngrid)
162            ix=int(xtn)
163            jy=int(ytn)
164            ddy=ytn-real(jy)
165            ddx=xtn-real(ix)
166          else
167            ix=int(xl)
168            jy=int(yl)
169            ddy=yl-real(jy)
170            ddx=xl-real(ix)
171          endif
172          ixp=ix+1
173          jyp=jy+1
174          rddx=1.-ddx
175          rddy=1.-ddy
176          p1=rddx*rddy
177          p2=ddx*rddy
178          p3=rddx*ddy
179          p4=ddx*ddy
180
181          if (ngrid.gt.0) then
182            oroh=oroh+p1*oron(ix ,jy ,ngrid) &
183                 + p2*oron(ixp,jy ,ngrid) &
184                 + p3*oron(ix ,jyp,ngrid) &
185                 + p4*oron(ixp,jyp,ngrid)
186          else
187            oroh=oroh+p1*oro(ix ,jy) &
188                 + p2*oro(ixp,jy) &
189                 + p3*oro(ix ,jyp) &
190                 + p4*oro(ixp,jyp)
191          endif
192        end do
193      end do
194
195  ! Divide by the number of samples taken
196  !**************************************
197
198      oroout(iix,jjy)=oroh/100.
199    end do
200  end do
201
202  ! if necessary allocate flux fields
203  if (iflux.eq.1) then
204    allocate(flux(6,0:numxgrid-1,0:numygrid-1,numzgrid, &
205         1:nspec,1:maxpointspec_act,1:nageclass),stat=stat)
206    if (stat.ne.0) write(*,*)'ERROR: could not allocate flux array '
207  endif
208
209  ! gridunc,griduncn        uncertainty of outputted concentrations
210  allocate(gridunc(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec, &
211       maxpointspec_act,nclassunc,maxageclass),stat=stat)
212    if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
213  if (ldirect.gt.0) then
214    allocate(wetgridunc(0:numxgrid-1,0:numygrid-1,maxspec, &
215       maxpointspec_act,nclassunc,maxageclass),stat=stat)
216    if (stat.ne.0) write(*,*)'ERROR: could not allocate wetgridunc'
217  allocate(drygridunc(0:numxgrid-1,0:numygrid-1,maxspec, &
218       maxpointspec_act,nclassunc,maxageclass),stat=stat)
219    if (stat.ne.0) write(*,*)'ERROR: could not allocate drygridunc'
220  endif
221
222! Extra field for totals at MPI root process
223  if (lroot.and.mpi_mode.gt.0) then
224
225#ifdef USE_MPIINPLACE
226#else
227    ! If MPI_IN_PLACE option is not used in mpi_mod.f90::mpif_tm_reduce_grid(),
228    ! then an aux array is needed for parallel grid reduction
229    allocate(gridunc0(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec, &
230         maxpointspec_act,nclassunc,maxageclass),stat=stat)
231    if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc0'
232#endif
233    if (ldirect.gt.0) then
234      allocate(wetgridunc0(0:numxgrid-1,0:numygrid-1,maxspec, &
235           maxpointspec_act,nclassunc,maxageclass),stat=stat)
236      if (stat.ne.0) write(*,*)'ERROR: could not allocate wetgridunc0'
237      allocate(drygridunc0(0:numxgrid-1,0:numygrid-1,maxspec, &
238           maxpointspec_act,nclassunc,maxageclass),stat=stat)
239      if (stat.ne.0) write(*,*)'ERROR: could not allocate drygridunc0'
240    endif
241! allocate a dummy to avoid compilator complaints
242  else if (.not.lroot.and.mpi_mode.gt.0) then
243    allocate(wetgridunc0(1,1,1,1,1,1),stat=stat)
244    allocate(drygridunc0(1,1,1,1,1,1),stat=stat)
245  end if
246
247  !write (*,*) 'Dimensions for fields', numxgrid,numygrid, &
248  !     maxspec,maxpointspec_act,nclassunc,maxageclass
249
250  if (lroot) then
251    write (*,*) 'Allocating fields for global output (x,y): ', numxgrid,numygrid
252    write (*,*) 'Allocating fields for nested output (x,y): ', numxgridn,numygridn
253  end if
254
255  ! allocate fields for concoutput with maximum dimension of outgrid
256  ! and outgrid_nest
257
258  allocate(gridsigma(0:max(numxgrid,numxgridn)-1, &
259       0:max(numygrid,numygridn)-1,numzgrid),stat=stat)
260    if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
261  allocate(grid(0:max(numxgrid,numxgridn)-1, &
262       0:max(numygrid,numygridn)-1,numzgrid),stat=stat)
263    if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
264  allocate(densityoutgrid(0:max(numxgrid,numxgridn)-1, &
265       0:max(numygrid,numygridn)-1,numzgrid),stat=stat)
266    if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
267
268  allocate(factor3d(0:max(numxgrid,numxgridn)-1, &
269       0:max(numygrid,numygridn)-1,numzgrid),stat=stat)
270    if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
271  allocate(sparse_dump_r(max(numxgrid,numxgridn)* &
272       max(numygrid,numygridn)*numzgrid),stat=stat)
273    if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
274
275   allocate(sparse_dump_u(max(numxgrid,numxgridn)* &
276       max(numygrid,numygridn)*numzgrid),stat=stat)
277        if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
278
279  allocate(sparse_dump_i(max(numxgrid,numxgridn)* &
280       max(numygrid,numygridn)*numzgrid),stat=stat)
281    if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
282
283  ! deposition fields are only allocated for forward runs
284  if (ldirect.gt.0) then
285     allocate(wetgridsigma(0:max(numxgrid,numxgridn)-1, &
286          0:max(numygrid,numygridn)-1),stat=stat)
287     if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
288     allocate(drygridsigma(0:max(numxgrid,numxgridn)-1, &
289          0:max(numygrid,numygridn)-1),stat=stat)
290     if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
291     allocate(wetgrid(0:max(numxgrid,numxgridn)-1, &
292          0:max(numygrid,numygridn)-1),stat=stat)
293     if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
294     allocate(drygrid(0:max(numxgrid,numxgridn)-1, &
295          0:max(numygrid,numygridn)-1),stat=stat)
296     if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
297  endif
298
299  ! Initial condition field
300
301  if (linit_cond.gt.0) then
302    allocate(init_cond(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec, &
303         maxpointspec_act),stat=stat)
304    if (stat.ne.0) write(*,*)'ERROR: could not allocate init_cond'
305  endif
306
307  !************************
308  ! Initialize output grids
309  !************************
310
311  do ks=1,nspec
312  do kp=1,maxpointspec_act
313    do i=1,numreceptor
314  ! Receptor points
315      creceptor(i,ks)=0.
316    end do
317    do nage=1,nageclass
318      do jy=0,numygrid-1
319        do ix=0,numxgrid-1
320          do l=1,nclassunc
321  ! Deposition fields
322            if (ldirect.gt.0) then
323            wetgridunc(ix,jy,ks,kp,l,nage)=0.
324            drygridunc(ix,jy,ks,kp,l,nage)=0.
325            endif
326            do kz=1,numzgrid
327              if (iflux.eq.1) then
328  ! Flux fields
329                 do i=1,5
330                   flux(i,ix,jy,kz,ks,kp,nage)=0.
331                 end do
332              endif
333  ! Initial condition field
334              if ((l.eq.1).and.(nage.eq.1).and.(linit_cond.gt.0)) &
335                   init_cond(ix,jy,kz,ks,kp)=0.
336  ! Concentration fields
337              gridunc(ix,jy,kz,ks,kp,l,nage)=0.
338            end do
339          end do
340        end do
341      end do
342    end do
343  end do
344  end do
345
346
347
348end subroutine outgrid_init
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG