source: trunk/src/assignland.f90 @ 28

Last change on this file since 28 was 4, checked in by mlanger, 11 years ago
File size: 9.3 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 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 :: xlandusep(lumaxx,lumaxy,numclass)
63  ! character*2 ck
64
65  do ix=1,lumaxx
66    do jy=1,lumaxy
67          do k=1,numclass
68            xlandusep(ix,jy,k)=0.
69          end do
70          sumperc=0.
71          do li=1,3
72            sumperc=sumperc+landinvent(ix,jy,li+3)
73          end do
74          do li=1,3
75            k=landinvent(ix,jy,li)
76          if (sumperc.gt.0) then
77            p=landinvent(ix,jy,li+3)/sumperc
78          else
79            p=0
80          endif
81  ! p has values between 0 and 1
82            xlandusep(ix,jy,k)=p
83          end do
84    end do
85  end do
86
87  ! do 13 k=1,11
88  ! write (ck,'(i2.2)') k
89  ! open(4,file='xlandusetest'//ck,form='formatted')
90  ! do 11 ix=1,lumaxx
91  !11       write (4,*) (xlandusep(ix,jy,k),jy=1,lumaxy)
92  !11       write (4,*) (landinvent(ix,jy,k),jy=1,lumaxy)
93  !13     close(4)
94
95  ! write (*,*) xlon0,ylat0,xlon0n(1),ylat0n(1),nxmin1,nymin1
96  ! write (*,*) dx, dy, dxout, dyout, ylat0, xlon0
97  nrefine=10
98  do ix=0,nxmin1
99    do jy=0,nymin1
100      do k=1,numclass
101        sumperc=0.
102        xlanduse(ix,jy,k)=0.
103      end do
104        do iix=1, nrefine
105          xlon=(ix+(iix-1)/real(nrefine))*dx+xlon0        ! longitude, should be between -180 and 179
106          if (xlon.ge.(xlon0lu+lumaxx*dxlu))  then
107               xlon=xlon-lumaxx*dxlu
108          endif
109          do jjy=1, nrefine
110           ylat=(jy+(jjy-1)/real(nrefine))*dy+ylat0       ! and lat. of each gridpoint
111           xi=int((xlon-xlon0lu)/dxlu)+1
112           yj=int((ylat-ylat0lu)/dxlu)+1
113           if (xi.gt.lumaxx) xi=xi-lumaxx
114           if (yj.gt.lumaxy) yj=yj-lumaxy
115           if (xi.lt.0) then
116              write (*,*) 'problem with landuseinv sampling: ', &
117                   xlon,xlon0lu,ix,iix,xlon0,dx,nxmax
118              stop
119           endif
120           do k=1,numclass
121              xlanduse(ix,jy,k)= &
122                   xlanduse(ix,jy,k)+xlandusep(int(xi),int(yj),k)
123             sumperc=sumperc+xlanduse(ix,jy,k)  ! just for the check if landuseinv. is available
124           end do
125          end do
126        end do
127          if (sumperc.gt.0) then                       ! detailed landuse available
128          sumperc=0.
129          do k=1,numclass
130              xlanduse(ix,jy,k)= &
131                   xlanduse(ix,jy,k)/real(nrefine*nrefine)
132            sumperc=sumperc+xlanduse(ix,jy,k)
133          end do
134  !cc the sum of all categories should be 1 ... 100 percent ... in order to get vdep right!
135          if (sumperc.lt.1-1E-5) then
136            do k=1,numclass
137              xlanduse(ix,jy,k)= &
138                   xlanduse(ix,jy,k)/sumperc
139            end do
140          endif
141          else
142            if (lsm(ix,jy).lt.0.1) then           ! over sea  -> ocean
143              xlanduse(ix,jy,3)=1.
144            else                                  ! over land -> rangeland
145              xlanduse(ix,jy,7)=1.
146            endif
147          endif
148
149
150    end do
151  end do
152
153  !***********************************
154  ! for test: write out xlanduse
155
156  ! open(4,file='landusetest',form='formatted')
157  ! do 56 k=1,13
158  ! do 55 ix=0,nxmin1
159  !55    write (4,*) (xlanduse(ix,jy,k),jy=0,nymin1)
160  !56    continue
161  ! close(4)
162  ! write (*,*) 'landuse written'
163  !stop
164  !  open(4,file='landseatest'//ck,form='formatted')
165  ! do 57 ix=0,nxmin1
166  !57       write (4,*) (lsm(ix,jy),jy=0,nymin1)
167  !  write (*,*) 'landseamask written'
168
169  !****************************************
170  ! Same as above, but for the nested grids
171  !****************************************
172
173  !************** TEST ********************
174  ! dyn(1)=dyn(1)/40
175  ! dxn(1)=dxn(1)/40
176  ! xlon0n(1)=1
177  ! ylat0n(1)=50
178  !************** TEST ********************
179
180  do l=1,numbnests
181    do ix=0,nxn(l)-1
182      do jy=0,nyn(l)-1
183        do k=1,numclass
184          sumperc=0.
185          xlandusen(ix,jy,k,l)=0.
186        end do
187          do iix=1, nrefine
188           xlon=(ix+(iix-1)/real(nrefine))*dxn(l)+xlon0n(l)
189           do jjy=1, nrefine
190             ylat=(jy+(jjy-1)/real(nrefine))*dyn(l)+ylat0n(l)
191             xi=int((xlon-xlon0lu)/dxlu)+1
192             yj=int((ylat-ylat0lu)/dxlu)+1
193             if (xi.gt.lumaxx) xi=xi-lumaxx
194             if (yj.gt.lumaxy) yj=yj-lumaxy
195             do k=1,numclass
196                xlandusen(ix,jy,k,l)=xlandusen(ix,jy,k,l)+ &
197                     xlandusep(int(xi),int(yj),k)
198               sumperc=sumperc+xlandusen(ix,jy,k,l)
199             end do
200           end do
201          end do
202          if (sumperc.gt.0) then                     ! detailed landuse available
203          sumperc=0.
204            do k=1,numclass
205               xlandusen(ix,jy,k,l)= &
206                    xlandusen(ix,jy,k,l)/real(nrefine*nrefine)
207              sumperc=sumperc+xlandusen(ix,jy,k,l)
208            end do
209  !cc the sum of all categories should be 1 ... 100 percent ... in order to get vdep right!
210            if (sumperc.lt.1-1E-5) then
211              do k=1,numclass
212                xlandusen(ix,jy,k,l)=xlandusen(ix,jy,k,l)/sumperc
213              end do
214            endif
215          else                                    ! check land/sea mask
216            if (lsmn(ix,jy,l).lt.0.1) then   ! over sea  -> ocean
217              xlandusen(ix,jy,3,l)=1.
218            else                        ! over land -> grasslands
219              xlandusen(ix,jy,7,l)=1.
220            endif
221          endif
222      end do
223    end do
224  end do
225
226  !***********************************
227  ! for test: write out xlanduse
228
229  ! do 66 k=1,11
230  ! write (ck,'(i2.2)') k
231  ! open(4,file='nlandusetest'//ck,form='formatted')
232  ! do 65 ix=0,nxn(1)-1
233  !65       write (4,*) (xlandusen(ix,jy,k,1),jy=0,nyn(1)-1)
234  !66      close(4)
235
236  ! write (*,*) 'landuse nested written'
237
238
239end subroutine assignland
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG