source: flexpart.git/src/outgrid_init_nest.f90 @ 8a65cb0

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bugunivie
Last change on this file since 8a65cb0 was 8a65cb0, checked in by Espen Sollum ATMOS <espen@…>, 9 years ago

Added code, makefile for dev branch

  • Property mode set to 100644
File size: 8.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_nest
23
24!*****************************************************************************
25!                                                                            *
26!  This routine calculates, for each grid cell of the output nest, the       *
27!  volume and the surface area.                                              *
28!                                                                            *
29!     Author: A. Stohl                                                       *
30!                                                                            *
31!    30 August 2004                                                          *
32!                                                                            *
33!*****************************************************************************
34!                                                                            *
35! Variables:                                                                 *
36!                                                                            *
37! arean              surface area of all output nest cells                   *
38! volumen            volumes of all output nest cells                        *
39!                                                                            *
40!*****************************************************************************
41
42  use unc_mod
43  use outg_mod
44  use par_mod
45  use com_mod
46
47  implicit none
48
49  integer :: ix,jy,kz,ks,kp,nage,l,iix,jjy,ixp,jyp,i1,j1,j,ngrid
50  integer :: stat
51  real :: ylat,gridarea,ylatp,ylatm,hzone,cosfactm,cosfactp
52  real :: xlon,xl,yl,ddx,ddy,rddx,rddy,p1,p2,p3,p4,xtn,ytn,oroh
53  real,parameter :: eps=nxmax/3.e5
54
55
56
57! gridunc,griduncn        uncertainty of outputted concentrations
58  allocate(griduncn(0:numxgridn-1,0:numygridn-1,numzgrid,maxspec, &
59       maxpointspec_act,nclassunc,maxageclass),stat=stat)
60  if (stat.ne.0) write(*,*)'ERROR:could not allocate nested gridunc'
61
62  if (ldirect.gt.0) then
63    allocate(wetgriduncn(0:numxgridn-1,0:numygridn-1,maxspec, &
64         maxpointspec_act,nclassunc,maxageclass),stat=stat)
65    if (stat.ne.0) write(*,*)'ERROR:could not allocate nested gridunc'
66    allocate(drygriduncn(0:numxgridn-1,0:numygridn-1,maxspec, &
67         maxpointspec_act,nclassunc,maxageclass),stat=stat)
68    if (stat.ne.0) write(*,*)'ERROR:could not allocate nested gridunc'
69  endif
70
71! Extra field for totals at MPI root process
72  if (lroot.and.mpi_mode.gt.0) then
73    ! allocate(griduncn0(0:numxgridn-1,0:numygridn-1,numzgrid,maxspec, &
74    !      maxpointspec_act,nclassunc,maxageclass),stat=stat)
75    ! if (stat.ne.0) write(*,*)'ERROR:could not allocate nested gridunc'
76
77    if (ldirect.gt.0) then
78      allocate(wetgriduncn0(0:numxgridn-1,0:numygridn-1,maxspec, &
79           maxpointspec_act,nclassunc,maxageclass),stat=stat)
80      if (stat.ne.0) write(*,*)'ERROR:could not allocate nested gridunc'
81      allocate(drygriduncn0(0:numxgridn-1,0:numygridn-1,maxspec, &
82           maxpointspec_act,nclassunc,maxageclass),stat=stat)
83      if (stat.ne.0) write(*,*)'ERROR:could not allocate nested gridunc'
84    endif
85! allocate a dummy to avoid compilator complaints
86  else if (.not.lroot.and.mpi_mode.gt.0) then
87    allocate(wetgriduncn0(1,1,1,1,1,1),stat=stat)
88    allocate(drygriduncn0(1,1,1,1,1,1),stat=stat)
89  end if
90
91! Compute surface area and volume of each grid cell: area, volume;
92! and the areas of the northward and eastward facing walls: areaeast, areanorth
93!***********************************************************************
94
95  do jy=0,numygridn-1
96    ylat=outlat0n+(real(jy)+0.5)*dyoutn
97    ylatp=ylat+0.5*dyoutn
98    ylatm=ylat-0.5*dyoutn
99    if ((ylatm.lt.0).and.(ylatp.gt.0.)) then
100      hzone=dyoutn*r_earth*pi180
101    else
102
103! Calculate area of grid cell with formula M=2*pi*R*h*dx/360,
104! see Netz, Formeln der Mathematik, 5. Auflage (1983), p.90
105!************************************************************
106
107      cosfactp=cos(ylatp*pi180)
108      cosfactm=cos(ylatm*pi180)
109      if (cosfactp.lt.cosfactm) then
110        hzone=sqrt(1-cosfactp**2)- &
111             sqrt(1-cosfactm**2)
112        hzone=hzone*r_earth
113      else
114        hzone=sqrt(1-cosfactm**2)- &
115             sqrt(1-cosfactp**2)
116        hzone=hzone*r_earth
117      endif
118    endif
119
120
121
122! Surface are of a grid cell at a latitude ylat
123!**********************************************
124
125    gridarea=2.*pi*r_earth*hzone*dxoutn/360.
126
127    do ix=0,numxgridn-1
128      arean(ix,jy)=gridarea
129
130! Volume = area x box height
131!***************************
132
133      volumen(ix,jy,1)=arean(ix,jy)*outheight(1)
134      do kz=2,numzgrid
135        volumen(ix,jy,kz)=arean(ix,jy)*(outheight(kz)-outheight(kz-1))
136      end do
137    end do
138  end do
139
140
141!**************************************************************************
142! Determine average height of model topography in nesteed output grid cells
143!**************************************************************************
144
145! Loop over all output grid cells
146!********************************
147
148  do jjy=0,numygridn-1
149    do iix=0,numxgridn-1
150      oroh=0.
151
152! Take 100 samples of the topography in every grid cell
153!******************************************************
154
155      do j1=1,10
156        ylat=outlat0n+(real(jjy)+real(j1)/10.-0.05)*dyoutn
157        yl=(ylat-ylat0)/dy
158        do i1=1,10
159          xlon=outlon0n+(real(iix)+real(i1)/10.-0.05)*dxoutn
160          xl=(xlon-xlon0)/dx
161
162! Determine the nest we are in
163!*****************************
164
165          ngrid=0
166          do j=numbnests,1,-1
167            if ((xl.gt.xln(j)+eps).and.(xl.lt.xrn(j)-eps).and. &
168                 (yl.gt.yln(j)+eps).and.(yl.lt.yrn(j)-eps)) then
169              ngrid=j
170              goto 43
171            endif
172          end do
17343        continue
174
175! Determine (nested) grid coordinates and auxiliary parameters used for interpolation
176!*****************************************************************************
177
178          if (ngrid.gt.0) then
179            xtn=(xl-xln(ngrid))*xresoln(ngrid)
180            ytn=(yl-yln(ngrid))*yresoln(ngrid)
181            ix=int(xtn)
182            jy=int(ytn)
183            ddy=ytn-real(jy)
184            ddx=xtn-real(ix)
185          else
186            ix=int(xl)
187            jy=int(yl)
188            ddy=yl-real(jy)
189            ddx=xl-real(ix)
190          endif
191          ixp=ix+1
192          jyp=jy+1
193          rddx=1.-ddx
194          rddy=1.-ddy
195          p1=rddx*rddy
196          p2=ddx*rddy
197          p3=rddx*ddy
198          p4=ddx*ddy
199
200          if (ngrid.gt.0) then
201            oroh=oroh+p1*oron(ix ,jy ,ngrid) &
202                 + p2*oron(ixp,jy ,ngrid) &
203                 + p3*oron(ix ,jyp,ngrid) &
204                 + p4*oron(ixp,jyp,ngrid)
205          else
206            oroh=oroh+p1*oro(ix ,jy) &
207                 + p2*oro(ixp,jy) &
208                 + p3*oro(ix ,jyp) &
209                 + p4*oro(ixp,jyp)
210          endif
211        end do
212      end do
213
214! Divide by the number of samples taken
215!**************************************
216
217      orooutn(iix,jjy)=oroh/100.
218    end do
219  end do
220
221
222
223!*******************************
224! Initialization of output grids
225!*******************************
226
227  do kp=1,maxpointspec_act
228    do ks=1,nspec
229      do nage=1,nageclass
230        do jy=0,numygridn-1
231          do ix=0,numxgridn-1
232            do l=1,nclassunc
233! Deposition fields
234              if (ldirect.gt.0) then
235                wetgriduncn(ix,jy,ks,kp,l,nage)=0.
236                drygriduncn(ix,jy,ks,kp,l,nage)=0.
237              endif
238! Concentration fields
239              do kz=1,numzgrid
240                griduncn(ix,jy,kz,ks,kp,l,nage)=0.
241              end do
242            end do
243          end do
244        end do
245      end do
246    end do
247  end do
248
249
250end subroutine outgrid_init_nest
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG