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

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

sources for flexwrf v3.1
File size: 10.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
24      subroutine outgrid_init_nest_reg
25!*******************************************************************************
26!                                                                              *
27!     Note:  This is the FLEXPART_WRF version of subroutine outgrid_init.      *
28!            The computational grid is the WRF x-y grid rather than lat-lon.   *
29!                                                                              *
30!  This routine calculates, for each grid cell of the output grid, the         *
31!  volume, the surface area, and the areas of the northward and eastward       *
32!  facing surfaces.                                                            *
33!                                                                              *
34!     Author: A. Stohl                                                         *
35!                                                                              *
36!     7 August 2002                                                            *
37!                                                                              *
38!    26 Oct 2005, R. Easter - changes in gridarea, areaeast, areanorth         *
39!                             associated with WRF horizontal grid.             *
40!     Dec 2005, R. Easter - changed names of "*lon0*" & "*lat0*" variables     *
41!     July 2012, J Brioude - modified for regular output grid.                 *
42!*******************************************************************************
43!                                                                              *
44! Variables:                                                                   *
45!                                                                              *
46! area               surface area of all output grid cells                     *
47! areaeast           eastward facing wall area of all output grid cells        *
48! areanorth          northward facing wall area of all output grid cells       *
49! volume             volumes of all output grid cells                          *
50!                                                                              *
51!*******************************************************************************
52
53  use unc_mod
54  use outg_mod
55  use par_mod
56  use com_mod
57
58!      include 'includepar'
59!      include 'includecom'
60     implicit none
61      integer :: ix,jy,kz,k,i,nage,l,iix,jjy,ixp,jyp,i1,j1,j,ngrid
62      real :: ylat,gridarea,ylatp,ylatm,hzone,cosfact,cosfactm,cosfactp
63      real :: ymet,xlon
64      real :: xmet,xl,yl,ddx,ddy,rddx,rddy,p1,p2,p3,p4,xtn,ytn,oroh
65  integer :: ks,kp,stat
66  real,parameter :: eps=nxmax/3.e5
67
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=outlat0n+(real(jy)+0.5)*dyoutln
77        ylatp=ylat+0.5*dyoutln
78        ylatm=ylat-0.5*dyoutln
79        if ((ylatm.lt.0).and.(ylatp.gt.0.)) then
80          hzone=dyoutln*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*dxoutln/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!       gridarea=dxoutn*dyoutn
108
109        do ix=0,numxgridn-1
110          arean(ix,jy)=gridarea
111
112! Volume = area x box height
113!***************************
114
115          volumen(ix,jy,1)=arean(ix,jy)*outheight(1)
116
117          do kz=2,numzgrid
118
119          volumen(ix,jy,kz)=arean(ix,jy)*(outheight(kz)-outheight(kz-1))
120      end do
121    end do
122  end do
123
124
125
126
127!******************************************************************
128! Determine average height of model topography in output grid cells
129!******************************************************************
130
131! Loop over all output grid cells
132!********************************
133
134      do jjy=0,numygridn-1
135        do iix=0,numxgridn-1
136          oroh=0.
137
138! Take 100 samples of the topography in every grid cell
139!******************************************************
140
141          do j1=1,10
142! for FLEXPART_WRF, x & y coords are in meters,
143! and the lon & lat variables below are in meters.
144            ylat=outlat0n+(float(jjy)+float(j1)/10.-0.05)*dyoutln !in degrees
145!           ymet=out_ym0n+(real(jjy)+real(j1)/10.-0.05)*dyoutn
146!           yl=(ymet-ymet0)/dy
147            do i1=1,10
148            xlon=outlon0n+(float(iix)+float(i1)/10.-0.05)*dxoutln !in degrees
149!              xmet=out_xm0n+(real(iix)+real(i1)/10.-0.05)*dxoutn
150        call ll_to_xymeter_wrf(xlon,ylat,xmet,ymet)
151           yl=(ymet-ymet0)/dy
152           xl=(xmet-xmet0)/dx
153!             xl=(xmet-xmet0)/dx
154
155! Determine the nest we are in
156!*****************************
157
158              ngrid=0
159              do j=numbnests,1,-1
160                if ((xl.gt.xln(j)).and.(xl.lt.xrn(j)).and. &
161                (yl.gt.yln(j)).and.(yl.lt.yrn(j))) then
162                  ngrid=j
163                  goto 43
164                endif
165          end do
16643            continue
167
168! Determine (nested) grid coordinates and auxiliary parameters used for interpolation
169!************************************************************************************
170
171              if (ngrid.gt.0) then
172                xtn=(xl-xln(ngrid))*xresoln(ngrid)
173                ytn=(yl-yln(ngrid))*yresoln(ngrid)
174                ix=int(xtn)
175                jy=int(ytn)
176                ddy=ytn-real(jy)
177                ddx=xtn-real(ix)
178              else
179                ix=int(xl)
180                jy=int(yl)
181                ddy=yl-real(jy)
182                ddx=xl-real(ix)
183              endif
184              ixp=ix+1
185              jyp=jy+1
186              rddx=1.-ddx
187              rddy=1.-ddy
188              p1=rddx*rddy
189              p2=ddx*rddy
190              p3=rddx*ddy
191              p4=ddx*ddy
192
193          if (ngrid.gt.0) then
194            oroh=oroh+p1*oron(ix ,jy ,ngrid) &
195                 + p2*oron(ixp,jy ,ngrid) &
196                 + p3*oron(ix ,jyp,ngrid) &
197                 + p4*oron(ixp,jyp,ngrid)
198          else
199            oroh=oroh+p1*oro(ix ,jy) &
200                 + p2*oro(ixp,jy) &
201                 + p3*oro(ix ,jyp) &
202                 + p4*oro(ixp,jyp)
203          endif
204        end do
205      end do
206
207! Divide by the number of samples taken
208!**************************************
209
210          orooutn(iix,jjy)=oroh/100.
211    end do
212  end do
213
214
215  ! gridunc,griduncn        uncertainty of outputted concentrations
216  allocate(griduncn(0:numxgridn-1,0:numygridn-1,numzgrid,maxspec, &
217       maxpointspec_act,nclassunc,maxageclass),stat=stat)
218  if (stat.ne.0) write(*,*)'ERROR:could not allocate nested gridunc'
219
220  if (ldirect.gt.0) then
221  allocate(wetgriduncn(0:numxgridn-1,0:numygridn-1,maxspec, &
222       maxpointspec_act,nclassunc,maxageclass),stat=stat)
223  if (stat.ne.0) write(*,*)'ERROR:could not allocate nested gridunc'
224  allocate(drygriduncn(0:numxgridn-1,0:numygridn-1,maxspec, &
225       maxpointspec_act,nclassunc,maxageclass),stat=stat)
226  allocate(drygriduncn2(0:numxgridn-1,0:numygridn-1,maxspec, &
227       maxpointspec_act,nclassunc,maxageclass),stat=stat)
228  if (stat.ne.0) write(*,*)'ERROR:could not allocate nested gridunc'
229  endif
230
231  !write (*,*) 'Dimensions for fields', numxgrid,numygrid, &
232  !     maxspec,maxpointspec_act,nclassunc,maxageclass
233
234  ! allocate fields for concoutput with maximum dimension of outgrid
235  ! and outgrid_nest
236  ! Initial condition field
237
238  !************************
239  ! Initialize output grids
240  !************************
241
242  do kp=1,maxpointspec_act
243  do ks=1,nspec
244    do nage=1,nageclass
245      do jy=0,numygridn-1
246        do ix=0,numxgridn-1
247          do l=1,nclassunc
248  ! Deposition fields
249            if (ldirect.gt.0) then
250              wetgriduncn(ix,jy,ks,kp,l,nage)=0.
251              drygriduncn(ix,jy,ks,kp,l,nage)=0.
252            endif
253  ! Concentration fields
254            do kz=1,numzgrid
255              griduncn(ix,jy,kz,ks,kp,l,nage)=0.
256            end do
257          end do
258        end do
259      end do
260    end do
261  end do
262  end do
263
264
265end subroutine outgrid_init_nest_reg
266
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG