assignland.f90 7.92 KB
Newer Older
Matthias Langer's avatar
 
Matthias Langer committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
subroutine assignland

  !*****************************************************************************
  !                                                                            *
  !     This routine assigns fractions of the 13 landuse classes to each ECMWF *
  !     grid point.                                                            *
  !     The landuse inventory of                                               *
  !                                                                            *
  ! Belward, A.S., Estes, J.E., and Kline, K.D., 1999,                         *
  ! The IGBP-DIS 1-Km Land-Cover Data Set DISCover:                            *
  ! A Project Overview: Photogrammetric Engineering and Remote Sensing ,       *
  ! v. 65, no. 9, p. 1013-1020                                                 *
  !                                                                            *
  !     if there are no data in the inventory                                  *
  !     the ECMWF land/sea mask is used to distinguish                         *
  !     between sea (-> ocean) and land (-> grasslands).                       *
  !                                                                            *
  !     Author: A. Stohl                                                       *
  !                                                                            *
  !     5 December 1996                                                        *
  !     8 February 1999 Additional use of nests, A. Stohl                      *
  !    29 December 2006 new landuse inventory, S. Eckhardt                     *
  !*****************************************************************************
  !                                                                            *
  ! Variables:                                                                 *
  ! xlanduse          fractions of numclass landuses for each model grid point *
  ! landinvent       landuse inventory (0.3 deg resolution)                    *
  !                                                                            *
  !*****************************************************************************

  use par_mod
  use com_mod

  implicit none

  integer :: ix,jy,k,l,li,nrefine,iix,jjy
  integer,parameter :: lumaxx=1200,lumaxy=600
  integer,parameter :: xlon0lu=-180,ylat0lu=-90
  real,parameter :: dxlu=0.3
  real :: xlon,ylat,sumperc,p,xi,yj
  real :: xlandusep(lumaxx,lumaxy,numclass)
  ! character*2 ck

  do ix=1,lumaxx
    do jy=1,lumaxy
          do k=1,numclass
            xlandusep(ix,jy,k)=0.
          end do
          sumperc=0.
          do li=1,3
            sumperc=sumperc+landinvent(ix,jy,li+3)
          end do
          do li=1,3
            k=landinvent(ix,jy,li)
          if (sumperc.gt.0) then
            p=landinvent(ix,jy,li+3)/sumperc
          else
            p=0
          endif
  ! p has values between 0 and 1
            xlandusep(ix,jy,k)=p
          end do
    end do
  end do

  ! do 13 k=1,11
  ! write (ck,'(i2.2)') k
  ! open(4,file='xlandusetest'//ck,form='formatted')
  ! do 11 ix=1,lumaxx
  !11       write (4,*) (xlandusep(ix,jy,k),jy=1,lumaxy)
  !11       write (4,*) (landinvent(ix,jy,k),jy=1,lumaxy)
  !13     close(4)

  ! write (*,*) xlon0,ylat0,xlon0n(1),ylat0n(1),nxmin1,nymin1
  ! write (*,*) dx, dy, dxout, dyout, ylat0, xlon0
  nrefine=10
  do ix=0,nxmin1
    do jy=0,nymin1
      do k=1,numclass
        sumperc=0.
        xlanduse(ix,jy,k)=0.
      end do
        do iix=1, nrefine
          xlon=(ix+(iix-1)/real(nrefine))*dx+xlon0        ! longitude, should be between -180 and 179
          if (xlon.ge.(xlon0lu+lumaxx*dxlu))  then
               xlon=xlon-lumaxx*dxlu
          endif
          do jjy=1, nrefine
           ylat=(jy+(jjy-1)/real(nrefine))*dy+ylat0       ! and lat. of each gridpoint
           xi=int((xlon-xlon0lu)/dxlu)+1
           yj=int((ylat-ylat0lu)/dxlu)+1
           if (xi.gt.lumaxx) xi=xi-lumaxx
           if (yj.gt.lumaxy) yj=yj-lumaxy
           if (xi.lt.0) then
              write (*,*) 'problem with landuseinv sampling: ', &
                   xlon,xlon0lu,ix,iix,xlon0,dx,nxmax
              stop
           endif
           do k=1,numclass
              xlanduse(ix,jy,k)= &
                   xlanduse(ix,jy,k)+xlandusep(int(xi),int(yj),k)
             sumperc=sumperc+xlanduse(ix,jy,k)  ! just for the check if landuseinv. is available
           end do
          end do
        end do
          if (sumperc.gt.0) then                       ! detailed landuse available
          sumperc=0.
          do k=1,numclass
              xlanduse(ix,jy,k)= &
                   xlanduse(ix,jy,k)/real(nrefine*nrefine)
            sumperc=sumperc+xlanduse(ix,jy,k)
          end do
  !cc the sum of all categories should be 1 ... 100 percent ... in order to get vdep right!
          if (sumperc.lt.1-1E-5) then
            do k=1,numclass
              xlanduse(ix,jy,k)= &
                   xlanduse(ix,jy,k)/sumperc
            end do
          endif
          else
            if (lsm(ix,jy).lt.0.1) then           ! over sea  -> ocean
              xlanduse(ix,jy,3)=1.
            else                                  ! over land -> rangeland
              xlanduse(ix,jy,7)=1.
            endif
          endif


    end do
  end do

  !***********************************
  ! for test: write out xlanduse

  ! open(4,file='landusetest',form='formatted')
  ! do 56 k=1,13
  ! do 55 ix=0,nxmin1
  !55    write (4,*) (xlanduse(ix,jy,k),jy=0,nymin1)
  !56    continue
  ! close(4)
  ! write (*,*) 'landuse written'
  !stop
  !  open(4,file='landseatest'//ck,form='formatted')
  ! do 57 ix=0,nxmin1
  !57       write (4,*) (lsm(ix,jy),jy=0,nymin1)
  !  write (*,*) 'landseamask written'

  !****************************************
  ! Same as above, but for the nested grids
  !****************************************

  !************** TEST ********************
  ! dyn(1)=dyn(1)/40
  ! dxn(1)=dxn(1)/40
  ! xlon0n(1)=1
  ! ylat0n(1)=50
  !************** TEST ********************

  do l=1,numbnests
    do ix=0,nxn(l)-1
      do jy=0,nyn(l)-1
        do k=1,numclass
          sumperc=0.
          xlandusen(ix,jy,k,l)=0.
        end do
          do iix=1, nrefine
           xlon=(ix+(iix-1)/real(nrefine))*dxn(l)+xlon0n(l)
           do jjy=1, nrefine
             ylat=(jy+(jjy-1)/real(nrefine))*dyn(l)+ylat0n(l)
             xi=int((xlon-xlon0lu)/dxlu)+1
             yj=int((ylat-ylat0lu)/dxlu)+1
             if (xi.gt.lumaxx) xi=xi-lumaxx
             if (yj.gt.lumaxy) yj=yj-lumaxy
             do k=1,numclass
                xlandusen(ix,jy,k,l)=xlandusen(ix,jy,k,l)+ &
                     xlandusep(int(xi),int(yj),k)
               sumperc=sumperc+xlandusen(ix,jy,k,l)
             end do
           end do
          end do
          if (sumperc.gt.0) then                     ! detailed landuse available
          sumperc=0.
            do k=1,numclass
               xlandusen(ix,jy,k,l)= &
                    xlandusen(ix,jy,k,l)/real(nrefine*nrefine)
              sumperc=sumperc+xlandusen(ix,jy,k,l)
            end do
  !cc the sum of all categories should be 1 ... 100 percent ... in order to get vdep right!
            if (sumperc.lt.1-1E-5) then
              do k=1,numclass
                xlandusen(ix,jy,k,l)=xlandusen(ix,jy,k,l)/sumperc
              end do
            endif
          else                                    ! check land/sea mask
            if (lsmn(ix,jy,l).lt.0.1) then   ! over sea  -> ocean
              xlandusen(ix,jy,3,l)=1.
            else                        ! over land -> grasslands
              xlandusen(ix,jy,7,l)=1.
            endif
          endif
      end do
    end do
  end do

  !***********************************
  ! for test: write out xlanduse

  ! do 66 k=1,11
  ! write (ck,'(i2.2)') k
  ! open(4,file='nlandusetest'//ck,form='formatted')
  ! do 65 ix=0,nxn(1)-1
  !65       write (4,*) (xlandusen(ix,jy,k,1),jy=0,nyn(1)-1)
  !66      close(4)

  ! write (*,*) 'landuse nested written'


end subroutine assignland