[e200b7a] | 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 :: 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 | |
---|
| 239 | end subroutine assignland |
---|