Context Navigation

assignland.f90

Last change on this file was 8, checked in by hasod, 11 years ago
Added parallel version of Flexpart91
File size: 9.5 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
22	subroutine assignland
23
24	!*****************************************************************************
25	! *
26	! This routine assigns fractions of the 13 landuse classes to each ECMWF *
27	! grid point. *
28	! The landuse inventory of *
29	! *
30	! Belward, A.S., Estes, J.E., and Kline, K.D., 1999, *
31	! The IGBP-DIS 1-Km Land-Cover Data Set DISCover: *
32	! A Project Overview: Photogrammetric Engineering and Remote Sensing , *
33	! v. 65, no. 9, p. 1013-1020 *
34	! *
35	! if there are no data in the inventory *
36	! the ECMWF land/sea mask is used to distinguish *
37	! between sea (-> ocean) and land (-> grasslands). *
38	! *
39	! Author: A. Stohl *
40	! *
41	! 5 December 1996 *
42	! 8 February 1999 Additional use of nests, A. Stohl *
43	! 29 December 2006 new landuse inventory, S. Eckhardt *
44	!*****************************************************************************
45	! *
46	! Variables: *
47	! xlanduse fractions of numclass landuses for each model grid point *
48	! landinvent landuse inventory (0.3 deg resolution) *
49	! *
50	!*****************************************************************************
51
52	use par_mod
53	use com_mod
54
55	implicit none
56
57	integer :: ix,jy,k,l,li,nrefine,iix,jjy
58	integer,parameter :: lumaxx=1200,lumaxy=600
59	integer,parameter :: xlon0lu=-180,ylat0lu=-90
60	real,parameter :: dxlu=0.3
61	real :: xlon,ylat,sumperc,p,xi,yj
62	real,allocatable,dimension(:,:,:) :: xlandusep
63	integer :: stat
64	! character*2 ck
65
66	allocate(xlandusep(lumaxx,lumaxy,numclass),stat=stat)
67	if (stat.ne.0) print*,'ERROR: could not allocate xlandusep'
68
69	do ix=1,lumaxx
70	do jy=1,lumaxy
71	do k=1,numclass
72	xlandusep(ix,jy,k)=0.
73	end do
74	sumperc=0.
75	do li=1,3
76	sumperc=sumperc+landinvent(ix,jy,li+3)
77	end do
78	do li=1,3
79	k=landinvent(ix,jy,li)
80	if (sumperc.gt.0) then
81	p=landinvent(ix,jy,li+3)/sumperc
82	else
83	p=0
84	endif
85	! p has values between 0 and 1
86	xlandusep(ix,jy,k)=p
87	end do
88	end do
89	end do
90
91	! do 13 k=1,11
92	! write (ck,'(i2.2)') k
93	! open(4,file='xlandusetest'//ck,form='formatted')
94	! do 11 ix=1,lumaxx
95	!11 write (4,*) (xlandusep(ix,jy,k),jy=1,lumaxy)
96	!11 write (4,*) (landinvent(ix,jy,k),jy=1,lumaxy)
97	!13 close(4)
98
99	! write (,) xlon0,ylat0,xlon0n(1),ylat0n(1),nxmin1,nymin1
100	! write (,) dx, dy, dxout, dyout, ylat0, xlon0
101	nrefine=10
102	do ix=0,nxmin1
103	do jy=0,nymin1
104	do k=1,numclass
105	sumperc=0.
106	xlanduse(ix,jy,k)=0.
107	end do
108	do iix=1, nrefine
109	xlon=(ix+(iix-1)/real(nrefine))*dx+xlon0 ! longitude, should be between -180 and 179
110	if (xlon.ge.(xlon0lu+lumaxx*dxlu)) then
111	xlon=xlon-lumaxx*dxlu
112	endif
113	do jjy=1, nrefine
114	ylat=(jy+(jjy-1)/real(nrefine))*dy+ylat0 ! and lat. of each gridpoint
115	xi=int((xlon-xlon0lu)/dxlu)+1
116	yj=int((ylat-ylat0lu)/dxlu)+1
117	if (xi.gt.lumaxx) xi=xi-lumaxx
118	if (yj.gt.lumaxy) yj=yj-lumaxy
119	if (xi.lt.0) then
120	write (,) 'problem with landuseinv sampling: ', &
121	xlon,xlon0lu,ix,iix,xlon0,dx,nxmax
122	stop
123	endif
124	do k=1,numclass
125	xlanduse(ix,jy,k)= &
126	xlanduse(ix,jy,k)+xlandusep(int(xi),int(yj),k)
127	sumperc=sumperc+xlanduse(ix,jy,k) ! just for the check if landuseinv. is available
128	end do
129	end do
130	end do
131	if (sumperc.gt.0) then ! detailed landuse available
132	sumperc=0.
133	do k=1,numclass
134	xlanduse(ix,jy,k)= &
135	xlanduse(ix,jy,k)/real(nrefine*nrefine)
136	sumperc=sumperc+xlanduse(ix,jy,k)
137	end do
138	!cc the sum of all categories should be 1 ... 100 percent ... in order to get vdep right!
139	if (sumperc.lt.1-1E-5) then
140	do k=1,numclass
141	xlanduse(ix,jy,k)= &
142	xlanduse(ix,jy,k)/sumperc
143	end do
144	endif
145	else
146	if (lsm(ix,jy).lt.0.1) then ! over sea -> ocean
147	xlanduse(ix,jy,3)=1.
148	else ! over land -> rangeland
149	xlanduse(ix,jy,7)=1.
150	endif
151	endif
152
153
154	end do
155	end do
156
157	!***********************************
158	! for test: write out xlanduse
159
160	! open(4,file='landusetest',form='formatted')
161	! do 56 k=1,13
162	! do 55 ix=0,nxmin1
163	!55 write (4,*) (xlanduse(ix,jy,k),jy=0,nymin1)
164	!56 continue
165	! close(4)
166	! write (,) 'landuse written'
167	!stop
168	! open(4,file='landseatest'//ck,form='formatted')
169	! do 57 ix=0,nxmin1
170	!57 write (4,*) (lsm(ix,jy),jy=0,nymin1)
171	! write (,) 'landseamask written'
172
173	!****************************************
174	! Same as above, but for the nested grids
175	!****************************************
176
177	!************ TEST ******************
178	! dyn(1)=dyn(1)/40
179	! dxn(1)=dxn(1)/40
180	! xlon0n(1)=1
181	! ylat0n(1)=50
182	!************ TEST ******************
183
184	do l=1,numbnests
185	do ix=0,nxn(l)-1
186	do jy=0,nyn(l)-1
187	do k=1,numclass
188	sumperc=0.
189	xlandusen(ix,jy,k,l)=0.
190	end do
191	do iix=1, nrefine
192	xlon=(ix+(iix-1)/real(nrefine))*dxn(l)+xlon0n(l)
193	do jjy=1, nrefine
194	ylat=(jy+(jjy-1)/real(nrefine))*dyn(l)+ylat0n(l)
195	xi=int((xlon-xlon0lu)/dxlu)+1
196	yj=int((ylat-ylat0lu)/dxlu)+1
197	if (xi.gt.lumaxx) xi=xi-lumaxx
198	if (yj.gt.lumaxy) yj=yj-lumaxy
199	do k=1,numclass
200	xlandusen(ix,jy,k,l)=xlandusen(ix,jy,k,l)+ &
201	xlandusep(int(xi),int(yj),k)
202	sumperc=sumperc+xlandusen(ix,jy,k,l)
203	end do
204	end do
205	end do
206	if (sumperc.gt.0) then ! detailed landuse available
207	sumperc=0.
208	do k=1,numclass
209	xlandusen(ix,jy,k,l)= &
210	xlandusen(ix,jy,k,l)/real(nrefine*nrefine)
211	sumperc=sumperc+xlandusen(ix,jy,k,l)
212	end do
213	!cc the sum of all categories should be 1 ... 100 percent ... in order to get vdep right!
214	if (sumperc.lt.1-1E-5) then
215	do k=1,numclass
216	xlandusen(ix,jy,k,l)=xlandusen(ix,jy,k,l)/sumperc
217	end do
218	endif
219	else ! check land/sea mask
220	if (lsmn(ix,jy,l).lt.0.1) then ! over sea -> ocean
221	xlandusen(ix,jy,3,l)=1.
222	else ! over land -> grasslands
223	xlandusen(ix,jy,7,l)=1.
224	endif
225	endif
226	end do
227	end do
228	end do
229
230	!***********************************
231	! for test: write out xlanduse
232
233	! do 66 k=1,11
234	! write (ck,'(i2.2)') k
235	! open(4,file='nlandusetest'//ck,form='formatted')
236	! do 65 ix=0,nxn(1)-1
237	!65 write (4,*) (xlandusen(ix,jy,k,1),jy=0,nyn(1)-1)
238	!66 close(4)
239
240	! write (,) 'landuse nested written'
241
242
243	end subroutine assignland

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format