source: flexpart.git/src/outgrid_init.f90

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bug
Last change on this file was 92fab65, checked in by Ignacio Pisso <ip@…>, 4 years ago

add SPDX-License-Identifier to all .f90 files

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