source: branches/jerome/src_flexwrf_v3.1/outgrid_init_nest_irreg.f90 @ 16

Last change on this file since 16 was 16, checked in by jebri, 10 years ago

sources for flexwrf v3.1

File size: 11.2 KB
Line 
1!***********************************************************************
2!* Copyright 2012,2013                                                *
3!* Jerome Brioude, Delia Arnold, Andreas Stohl, Wayne Angevine,       *
4!* John Burkhart, Massimo Cassiani, Adam Dingwell, Richard C Easter, Sabine Eckhardt,*
5!* Stephanie Evan, Jerome D Fast, Don Morton, Ignacio Pisso,          *
6!* Petra Seibert, Gerard Wotawa, Caroline Forster, Harald Sodemann,   *
7!*                                                                     *
8!* This file is part of FLEXPART WRF                                   *
9!*                                                                     *
10!* FLEXPART is free software: you can redistribute it and/or modify    *
11!* it under the terms of the GNU General Public License as published by*
12!* the Free Software Foundation, either version 3 of the License, or   *
13!* (at your option) any later version.                                 *
14!*                                                                     *
15!* FLEXPART is distributed in the hope that it will be useful,         *
16!* but WITHOUT ANY WARRANTY; without even the implied warranty of      *
17!* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
18!* GNU General Public License for more details.                        *
19!*                                                                     *
20!* You should have received a copy of the GNU General Public License   *
21!* along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
22!***********************************************************************
23      subroutine outgrid_init_nest_irreg
24!*******************************************************************************
25!                                                                              *
26!     Note:  This is the FLEXPART_WRF version of subroutine outgrid_init.      *
27!            The computational grid is the WRF x-y grid rather than lat-lon.   *
28!                                                                              *
29!  This routine calculates, for each grid cell of the output grid, the         *
30!  volume, the surface area, and the areas of the northward and eastward       *
31!  facing surfaces.                                                            *
32!                                                                              *
33!     Author: A. Stohl                                                         *
34!                                                                              *
35!     7 August 2002                                                            *
36!                                                                              *
37!    26 Oct 2005, R. Easter - changes in gridarea, areaeast, areanorth         *
38!                             associated with WRF horizontal grid.             *
39!     Dec 2005, R. Easter - changed names of "*lon0*" & "*lat0*" variables     *
40!                                                                              *
41!*******************************************************************************
42!                                                                              *
43! Variables:                                                                   *
44!                                                                              *
45! area               surface area of all output grid cells                     *
46! areaeast           eastward facing wall area of all output grid cells        *
47! areanorth          northward facing wall area of all output grid cells       *
48! volume             volumes of all output grid cells                          *
49!                                                                              *
50!*******************************************************************************
51
52  use unc_mod
53  use outg_mod
54  use par_mod
55  use com_mod
56
57!      include 'includepar'
58!      include 'includecom'
59     implicit none
60      integer :: ix,jy,kz,k,i,nage,l,iix,jjy,ixp,jyp,i1,j1,j,ngrid
61!     real ylat,gridarea,ylatp,ylatm,hzone,cosfact,cosfactm,cosfactp
62      real :: ymet,gridarea,m1,m2,xl1,xl2,yl1,yl2,tmpx,tmpy
63      real :: xmet,xl,yl,ddx,ddy,rddx,rddy,p1,p2,p3,p4,xtn,ytn,oroh
64  integer :: ks,kp,stat
65  real,parameter :: eps=nxmax/3.e5
66      real :: lon2(4),lat2(4)
67      real ( kind = 8 ) :: sphere01_polygon_area,haversine,area1
68
69
70! Compute surface area and volume of each grid cell: area, volume;
71! and the areas of the northward and eastward facing walls: areaeast, areanorth
72!***********************************************************************
73
74      do jy=0,numygridn-1
75
76!        ylat=outlat0+(real(jy)+0.5)*dyout
77!        ylatp=ylat+0.5*dyout
78!        ylatm=ylat-0.5*dyout
79!        if ((ylatm.lt.0).and.(ylatp.gt.0.)) then
80!          hzone=dyout*r_earth*pi180
81!        else
82!
83!C Calculate area of grid cell with formula M=2*pi*R*h*dx/360,
84!C see Netz, Formeln der Mathematik, 5. Auflage (1983), p.90
85!*************************************************************
86!
87!          cosfact=cos(ylat*pi180)*r_earth
88!          cosfactp=cos(ylatp*pi180)*r_earth
89!          cosfactm=cos(ylatm*pi180)*r_earth
90!          if (cosfactp.lt.cosfactm) then
91!            hzone=sqrt(r_earth**2-cosfactp**2)-
92!     +      sqrt(r_earth**2-cosfactm**2)
93!          else
94!            hzone=sqrt(r_earth**2-cosfactm**2)-
95!     +      sqrt(r_earth**2-cosfactp**2)
96!          endif
97!        endif
98!
99!C Surface are of a grid cell at a latitude ylat
100!***********************************************
101!
102!        gridarea=2.*pi*r_earth*hzone*dxout/360.
103
104! for FLEXPART_WRF, dx & dy are in meters, and no cos(lat) is needed
105! ??? maybe should incorporate map factor here,
106!     and also for areaeast & areanorth ???
107
108!       gridarea=dxoutn*dyoutn
109
110        do ix=0,numxgridn-1
111!        xl1=(real(ix)*dxoutn+out_xm0n)/dx
112!        yl1=(real(jy)*dyoutn+out_ym0n)/dy
113!        xl2=(real(ix+1)*dxoutn+out_xm0n)/dx
114!        yl2=(real(jy+1)*dyoutn+out_ym0n)/dy
115!!      xr=out_xm0+real(numxgrid)*dxout
116!        m1=0.5*(m_x(int(xl1),int(yl1),1)+m_x(int(xl2),int(yl1),1))
117!        m2=0.5*(m_y(int(xl1),int(yl1),1)+m_x(int(xl1),int(yl2),1))
118!      arean(ix,jy)=dxoutn*m1*dyoutn*m2
119
120! A more precise method
121          tmpx=out_xm0n+(float(ix))*dxoutn
122          tmpy=out_ym0n+(float(jy))*dyoutn
123          call xymeter_to_ll_wrf_out(tmpx,tmpy,lon2(1),lat2(1))
124          tmpx=out_xm0n+(float(ix+1))*dxoutn
125          tmpy=out_ym0n+(float(jy))*dyoutn
126          call xymeter_to_ll_wrf_out(tmpx,tmpy,lon2(2),lat2(2))
127          tmpx=out_xm0n+(float(ix+1))*dxoutn
128          tmpy=out_ym0n+(float(jy+1))*dyoutn
129          call xymeter_to_ll_wrf_out(tmpx,tmpy,lon2(3),lat2(3))
130          tmpx=out_xm0n+(float(ix))*dxoutn
131          tmpy=out_ym0n+(float(jy+1))*dyoutn
132          call xymeter_to_ll_wrf_out(tmpx,tmpy,lon2(4),lat2(4))
133        area1=sphere01_polygon_area ( 4, real(lat2,kind=8), real(lon2,kind=8) )
134        arean(ix,jy)=real(area1)*6370000.*6370000./coefdx/coefdx
135
136!         arean(ix,jy)=gridarea
137
138! Volume = area x box height
139!***************************
140
141          volumen(ix,jy,1)=arean(ix,jy)*outheight(1)
142
143          do kz=2,numzgrid
144
145          volumen(ix,jy,kz)=arean(ix,jy)*(outheight(kz)-outheight(kz-1))
146      end do
147    end do
148  end do
149
150
151
152
153!******************************************************************
154! Determine average height of model topography in output grid cells
155!******************************************************************
156
157! Loop over all output grid cells
158!********************************
159
160      do jjy=0,numygridn-1
161        do iix=0,numxgridn-1
162          oroh=0.
163
164! Take 100 samples of the topography in every grid cell
165!******************************************************
166
167          do j1=1,10
168! for FLEXPART_WRF, x & y coords are in meters,
169! and the lon & lat variables below are in meters.
170            ymet=out_ym0n+(real(jjy)+real(j1)/10.-0.05)*dyoutn
171            yl=(ymet-ymet0)/dy
172            do i1=1,10
173              xmet=out_xm0n+(real(iix)+real(i1)/10.-0.05)*dxoutn
174              xl=(xmet-xmet0)/dx
175
176! Determine the nest we are in
177!*****************************
178
179              ngrid=0
180              do j=numbnests,1,-1
181                if ((xl.gt.xln(j)).and.(xl.lt.xrn(j)).and. &
182                (yl.gt.yln(j)).and.(yl.lt.yrn(j))) then
183                  ngrid=j
184                  goto 43
185                endif
186          end do
18743            continue
188
189! Determine (nested) grid coordinates and auxiliary parameters used for interpolation
190!************************************************************************************
191
192              if (ngrid.gt.0) then
193                xtn=(xl-xln(ngrid))*xresoln(ngrid)
194                ytn=(yl-yln(ngrid))*yresoln(ngrid)
195                ix=int(xtn)
196                jy=int(ytn)
197                ddy=ytn-real(jy)
198                ddx=xtn-real(ix)
199              else
200                ix=int(xl)
201                jy=int(yl)
202                ddy=yl-real(jy)
203                ddx=xl-real(ix)
204              endif
205              ixp=ix+1
206              jyp=jy+1
207              rddx=1.-ddx
208              rddy=1.-ddy
209              p1=rddx*rddy
210              p2=ddx*rddy
211              p3=rddx*ddy
212              p4=ddx*ddy
213
214          if (ngrid.gt.0) then
215            oroh=oroh+p1*oron(ix ,jy ,ngrid) &
216                 + p2*oron(ixp,jy ,ngrid) &
217                 + p3*oron(ix ,jyp,ngrid) &
218                 + p4*oron(ixp,jyp,ngrid)
219          else
220            oroh=oroh+p1*oro(ix ,jy) &
221                 + p2*oro(ixp,jy) &
222                 + p3*oro(ix ,jyp) &
223                 + p4*oro(ixp,jyp)
224          endif
225        end do
226      end do
227
228! Divide by the number of samples taken
229!**************************************
230
231          orooutn(iix,jjy)=oroh/100.
232    end do
233  end do
234
235
236  ! gridunc,griduncn        uncertainty of outputted concentrations
237  allocate(griduncn(0:numxgridn-1,0:numygridn-1,numzgrid,maxspec, &
238       maxpointspec_act,nclassunc,maxageclass),stat=stat)
239  if (stat.ne.0) write(*,*)'ERROR:could not allocate nested gridunc'
240
241  if (ldirect.gt.0) then
242  allocate(wetgriduncn(0:numxgridn-1,0:numygridn-1,maxspec, &
243       maxpointspec_act,nclassunc,maxageclass),stat=stat)
244  if (stat.ne.0) write(*,*)'ERROR:could not allocate nested gridunc'
245  allocate(drygriduncn(0:numxgridn-1,0:numygridn-1,maxspec, &
246       maxpointspec_act,nclassunc,maxageclass),stat=stat)
247  allocate(drygriduncn2(0:numxgridn-1,0:numygridn-1,maxspec, &
248       maxpointspec_act,nclassunc,maxageclass),stat=stat)
249  if (stat.ne.0) write(*,*)'ERROR:could not allocate nested gridunc'
250  endif
251
252  !write (*,*) 'Dimensions for fields', numxgrid,numygrid, &
253  !     maxspec,maxpointspec_act,nclassunc,maxageclass
254
255  ! allocate fields for concoutput with maximum dimension of outgrid
256  ! and outgrid_nest
257  ! Initial condition field
258
259  !************************
260  ! Initialize output grids
261  !************************
262
263  do kp=1,maxpointspec_act
264  do ks=1,nspec
265    do nage=1,nageclass
266      do jy=0,numygridn-1
267        do ix=0,numxgridn-1
268          do l=1,nclassunc
269  ! Deposition fields
270            if (ldirect.gt.0) then
271              wetgriduncn(ix,jy,ks,kp,l,nage)=0.
272              drygriduncn(ix,jy,ks,kp,l,nage)=0.
273            endif
274  ! Concentration fields
275            do kz=1,numzgrid
276              griduncn(ix,jy,kz,ks,kp,l,nage)=0.
277            end do
278          end do
279        end do
280      end do
281    end do
282  end do
283  end do
284
285
286end subroutine outgrid_init_nest_irreg
287
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG