source: branches/jerome/src_flexwrf_v3.1/assignland.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: 9.8 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 assignland
25!*******************************************************************************
26!                                                                              *
27!     Note:  This is the FLEXPART_WRF version of subroutine assignland.        *
28!            The computational grid is the WRF x-y grid rather than lat-lon.   *
29!                                                                              *
30!     This routine assigns fractions of the 8 landuse classes to each ECMWF    *
31!     grid point.                                                              *
32!     The landuse inventory of                                                 *
33!                                                                              *
34!     van de Velde R.J., Faber W.S., van Katwijk V.F., Kuylenstierna J.C.I.,   *
35!     Scholten H.J., Thewessen T.J.M., Verspuij M., Zevenbergen M. (1994):     *
36!     The Preparation of a European Land Use Database. National Institute of   *
37!     Public Health and Environmental Protection, Report nr 712401001,         *
38!     Bilthoven, The Netherlands.                                              *
39!                                                                              *
40!     is used to create a detailed landuse inventory for Europe.               *
41!                                                                              *
42!     Outside of Europe, the ECMWF land/sea mask is used to distinguish        *
43!     between sea (-> ocean) and land (-> grasslands).                         *
44!                                                                              *
45!     Author: A. Stohl                                                         *
46!                                                                              *
47!     5 December 1996                                                          *
48!     8 February 1999 Additional use of nests, A. Stohl                        *
49!                                                                              *
50!    14 October  2005 R. Easter -- modified for WRF.                           *
51!                     The landuse inventory is not used at all.                *
52!                     The land/sea mask is used everywhere.                    *
53!    Feb 2012 J. Brioude: modified to fortran 90                               *
54!*******************************************************************************
55!                                                                              *
56! Variables:                                                                   *
57! xlanduse          fractions of numclass landuses for each model grid point   *
58! xlandinvent       landuse inventory (1/6 deg resolution)                     *
59!                                                                              *
60!*******************************************************************************
61
62  use par_mod
63  use com_mod
64
65  implicit none
66  integer :: ix,jy,k,l,li,nrefine,iix,jjy,n,nxn2,nyn2
67  integer,parameter :: lumaxx=1200,lumaxy=600
68  real,parameter :: xlon0lu=-180.,ylat0lu=-90.
69  real,parameter :: dxlu=0.3
70  real :: xlon,ylat,sumperc,p,xi,yj
71  real :: xlandusep(lumaxx,lumaxy,numclass)
72!  real :: xlanduse(0:nxmax-1,0:nymax-1,numclass)
73!  real :: xlandusen(0:nxmaxn-1,0:nymaxn-1,numclass,maxnests)
74
75!      include 'includepar'
76!      include 'includecom'
77!
78!      integer ix,jy,i,j,k,n,l
79!      real x,y,xlon,ylat
80
81  do ix=1,lumaxx
82    do jy=1,lumaxy
83          do k=1,numclass
84            xlandusep(ix,jy,k)=0.
85          end do
86          sumperc=0.
87          do li=1,3
88            sumperc=sumperc+landinvent(ix,jy,li+3)
89          end do
90          do li=1,3
91            k=landinvent(ix,jy,li)
92          if (sumperc.gt.0) then
93            p=landinvent(ix,jy,li+3)/sumperc
94          else
95            p=0
96          endif
97  ! p has values between 0 and 1
98            xlandusep(ix,jy,k)=p
99          end do
100    end do
101  end do
102
103  nrefine=10
104
105      do ix=0,nxmin1
106!       xlon=real(ix)*dx+xlon0         ! long.
107        do jy=0,nymin1
108!         ylat=real(jy)*dy+ylat0       ! and lat. of each gridpoint
109! FLEXPART_WRF - use this routine to get lat,lon
110          do k=1,numclass
111        sumperc=0.
112          xlanduse(ix,jy,k)=0.
113         enddo
114
115!        do iix=1, nrefine
116!          do jjy=1, nrefine
117        do iix=1, 1
118          do jjy=1, 1
119
120!          xlon=(ix+(iix-1)/real(nrefine))*dx+xlon0        ! longitude, should be between -180 and 179
121!           ylat=(jy+(jjy-1)/real(nrefine))*dy+ylat0       ! and lat. of each gridpoint
122
123          call xyindex_to_ll_wrf( 0, real(ix), real(jy), xlon, ylat )
124
125          if (xlon.ge.(xlon0lu+lumaxx*dxlu))  then
126               xlon=xlon-lumaxx*dxlu
127          endif
128           xi=int((xlon-xlon0lu)/dxlu)+1
129           yj=int((ylat-ylat0lu)/dxlu)+1
130           if (xi.gt.lumaxx) xi=xi-lumaxx
131           if (yj.gt.lumaxy) yj=yj-lumaxy
132           if (xi.lt.0) then
133              write (*,*) 'problem with landuseinv sampling: ', &
134                   xlon,xlon0lu,ix,iix,xlon0,dx,nxmax
135              stop
136           endif
137           do k=1,numclass
138              xlanduse(ix,jy,k)= &
139                   xlanduse(ix,jy,k)+xlandusep(int(xi),int(yj),k)
140             sumperc=sumperc+xlanduse(ix,jy,k)  ! just for the check if landuseinv. is available
141           end do
142          end do
143        end do
144          if (sumperc.gt.0) then                       ! detailed landuse available
145          sumperc=0.
146          do k=1,numclass
147              xlanduse(ix,jy,k)= &
148                   xlanduse(ix,jy,k)/real(nrefine*nrefine)
149            sumperc=sumperc+xlanduse(ix,jy,k)
150          end do
151  !cc the sum of all categories should be 1 ... 100 percent ... in order to get vdep right!
152          if (sumperc.lt.1-1E-5) then
153            do k=1,numclass
154              xlanduse(ix,jy,k)= &
155                   xlanduse(ix,jy,k)/sumperc
156            end do
157          endif
158          else
159            if (lsm(ix,jy).lt.0.1) then           ! over sea  -> ocean
160              xlanduse(ix,jy,3)=1.
161            else                                  ! over land -> rangeland
162              xlanduse(ix,jy,7)=1.
163            endif
164          endif
165
166
167    end do
168  end do
169
170
171!   open(4,file='landusetest',form='formatted')
172!   do 56 k=1,13
173 !  do 55 ix=0,nxmin1
174!  55    write (4,*) (xlanduse(ix,jy,k),jy=0,nymin1)
175!  56    continue
176!   close(4)
177!   write (*,*) 'landuse written'
178 ! stop
179 !   open(4,file='landseatest'//ck,form='formatted')
180 !  do 57 ix=0,nxmin1
181!  57       write (4,*) (lsm(ix,jy),jy=0,nymin1)
182!    write (*,*) 'landseamask written'
183
184
185      do l=1,numbnests
186       do jy=0,nyn(l)-1
187        do ix=0,nxn(l)-1
188
189          do k=1,numclass
190           xlandusen(ix,jy,k,l)=0.0
191          sumperc=0.
192           enddo
193
194        do iix=1, 1
195          do jjy=1, 1
196
197!          xlon=(ix+(iix-1)/real(nrefine))*dx+xlon0        ! longitude, should be between -180 and 179
198!           ylat=(jy+(jjy-1)/real(nrefine))*dy+ylat0       ! and lat. of each gridpoint
199
200          call xyindex_to_ll_wrf( l, real(ix), real(jy), xlon, ylat )
201
202             xi=int((xlon-xlon0lu)/dxlu)+1
203             yj=int((ylat-ylat0lu)/dxlu)+1
204             if (xi.gt.lumaxx) xi=xi-lumaxx
205             if (yj.gt.lumaxy) yj=yj-lumaxy
206             do k=1,numclass
207                xlandusen(ix,jy,k,l)=xlandusen(ix,jy,k,l)+ &
208                     xlandusep(int(xi),int(yj),k)
209               sumperc=sumperc+xlandusen(ix,jy,k,l)
210             end do
211           end do
212          end do
213          if (sumperc.gt.0) then                     ! detailed landuse available
214          sumperc=0.
215            do k=1,numclass
216               xlandusen(ix,jy,k,l)= &
217                    xlandusen(ix,jy,k,l)/real(nrefine*nrefine)
218              sumperc=sumperc+xlandusen(ix,jy,k,l)
219            end do
220  !cc the sum of all categories should be 1 ... 100 percent ... in order to get vdep right!
221            if (sumperc.lt.1-1E-5) then
222              do k=1,numclass
223                xlandusen(ix,jy,k,l)=xlandusen(ix,jy,k,l)/sumperc
224              end do
225            endif
226          else                                    ! check land/sea mask
227            if (lsmn(ix,jy,l).lt.0.1) then   ! over sea  -> ocean
228              xlandusen(ix,jy,3,l)=1.
229            else                        ! over land -> grasslands
230              xlandusen(ix,jy,7,l)=1.
231            endif
232          endif
233      end do
234    end do
235  end do
236
237
238      end subroutine assignland
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG