source: flexpart.git/src/assignland.f90 @ 3481cc1

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bug
Last change on this file since 3481cc1 was 3481cc1, checked in by Ignacio Pisso <ip@…>, 4 years ago

move license from headers to a different file

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