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