Changeset 37 for branches/petra/src
- Timestamp:
- Apr 22, 2015, 3:42:35 PM (9 years ago)
- Location:
- branches/petra/src
- Files:
-
- 1 added
- 17 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/petra/src/FLEXPART.f90
r36 r37 66 66 67 67 ! FLEXPART version string 68 flexversion='Version 9.2_beta.r3 4 (2015-02-12)'68 flexversion='Version 9.2_beta.r3? (2015-02-25)' 69 69 verbosity=0 70 70 -
branches/petra/src/com_mod.f90
r36 r37 32 32 ! Modifications: 15 August 2013 IP, 33 33 ! 2/2015 PS, add incremental deposition switch 34 ! 2/2015 PS, change variables for cloud diagnostics 34 35 ! * 35 36 !******************************************************************************* … … 365 366 real :: qvh(0:nxmax-1,0:nymax-1,nuvzmax,2) 366 367 real :: pplev(0:nxmax-1,0:nymax-1,nuvzmax,2) 367 integer(kind=1) :: clouds(0:nxmax-1,0:nymax-1,nzmax,2)368 integer :: cloudsh(0:nxmax-1,0:nymax-1,2)369 370 371 368 ! uu,vv,ww [m/2] wind components in x,y and z direction 372 369 ! uupol,vvpol [m/s] wind components in polar stereographic projection … … 377 374 ! drhodz [kg/m2] vertical air density gradient 378 375 ! tth,qvh tth,qvh on original eta levels 379 ! clouds: no cloud, no precipitation 0380 ! cloud, no precipitation 1381 ! rainout conv/lsp dominated 2/3382 ! washout conv/lsp dominated 4/5383 376 384 377 ! pplev for the GFS version … … 406 399 real :: oli(0:nxmax-1,0:nymax-1,1,2) 407 400 real :: diffk(0:nxmax-1,0:nymax-1,1,2) 401 integer, dimension (0:nxmax-1,0:nymax-1,2) :: icloudbot, icloudthck 408 402 409 403 ! ps surface pressure … … 426 420 ! oli [m] inverse Obukhov length (1/L) 427 421 ! diffk [m2/s] diffusion coefficient at reference height 422 ! icloudbot (m) cloud bottom height 423 ! icloudthck (m) cloud thickness 428 424 429 425 … … 484 480 real :: qvn(0:nxmaxn-1,0:nymaxn-1,nzmax,2,maxnests) 485 481 real :: pvn(0:nxmaxn-1,0:nymaxn-1,nzmax,2,maxnests) 486 integer(kind=1) :: cloudsn(0:nxmaxn-1,0:nymaxn-1,0:nzmax,2,maxnests)487 integer :: cloudsnh(0:nxmaxn-1,0:nymaxn-1,2,maxnests)488 482 real :: rhon(0:nxmaxn-1,0:nymaxn-1,nzmax,2,maxnests) 489 483 real :: drhodzn(0:nxmaxn-1,0:nymaxn-1,nzmax,2,maxnests) … … 514 508 real :: diffkn(0:nxmaxn-1,0:nymaxn-1,1,2,maxnests) 515 509 real :: vdepn(0:nxmaxn-1,0:nymaxn-1,maxspec,2,maxnests) 510 integer, dimension (0:nxmax-1,0:nymax-1,2,maxnests) :: icloudbotn, icloudthckn 516 511 517 512 -
branches/petra/src/getfields.f90
r4 r37 41 41 ! Variables tth,qvh,tthn,qvhn (on eta coordinates) in common block. 42 42 ! Function of nstop extended. 43 ! 44 ! 3/2015, PS: put in verbosity output 45 ! 43 46 !***************************************************************************** 44 47 ! * … … 83 86 integer :: indmin = 1 84 87 88 call verboutput(verbosity,1,' getfields') 85 89 86 90 ! Check, if wind fields are available for the current time step … … 125 129 do indj=indmin,numbwf-1 126 130 if (ldirect*wftime(indj+1).gt.ldirect*itime) then 131 call verboutput(verbosity,1,' call readwind') 127 132 call readwind(indj+1,memind(2),uuh,vvh,wwh) 133 call verboutput(verbosity,1,' call readwind_nests') 128 134 call readwind_nests(indj+1,memind(2),uuhn,vvhn,wwhn) 135 call verboutput(verbosity,1,' call calcpar') 129 136 call calcpar(memind(2),uuh,vvh,pvh) 137 call verboutput(verbosity,1,' call calcpar_nests') 130 138 call calcpar_nests(memind(2),uuhn,vvhn,pvhn) 139 call verboutput(verbosity,1,' call verttransform') 131 140 call verttransform(memind(2),uuh,vvh,wwh,pvh) 141 call verboutput(verbosity,1,' call verttransform_nests') 132 142 call verttransform_nests(memind(2),uuhn,vvhn,wwhn,pvhn) 133 143 memtime(2)=wftime(indj+1) … … 148 158 (ldirect*wftime(indj+1).gt.ldirect*itime)) then 149 159 memind(1)=1 160 call verboutput(verbosity,1,' call readwind') 150 161 call readwind(indj,memind(1),uuh,vvh,wwh) 162 call verboutput(verbosity,1,' call readwind_nests') 151 163 call readwind_nests(indj,memind(1),uuhn,vvhn,wwhn) 164 call verboutput(verbosity,1,' call calcpar') 152 165 call calcpar(memind(1),uuh,vvh,pvh) 166 call verboutput(verbosity,1,' call calcpar_nests') 153 167 call calcpar_nests(memind(1),uuhn,vvhn,pvhn) 168 call verboutput(verbosity,1,' call verttransform') 154 169 call verttransform(memind(1),uuh,vvh,wwh,pvh) 170 call verboutput(verbosity,1,' call verttransform_nests') 155 171 call verttransform_nests(memind(1),uuhn,vvhn,wwhn,pvhn) 156 172 memtime(1)=wftime(indj) 157 173 memind(2)=2 174 call verboutput(verbosity,1,' call readwind') 158 175 call readwind(indj+1,memind(2),uuh,vvh,wwh) 176 call verboutput(verbosity,1,' call readwind_nests') 159 177 call readwind_nests(indj+1,memind(2),uuhn,vvhn,wwhn) 178 call verboutput(verbosity,1,' call calcpar') 160 179 call calcpar(memind(2),uuh,vvh,pvh) 180 call verboutput(verbosity,1,' call calcpar_nests') 161 181 call calcpar_nests(memind(2),uuhn,vvhn,pvhn) 182 call verboutput(verbosity,1,' call verttransform') 162 183 call verttransform(memind(2),uuh,vvh,wwh,pvh) 184 call verboutput(verbosity,1,' call verttransform_nests') 163 185 call verttransform_nests(memind(2),uuhn,vvhn,wwhn,pvhn) 164 186 memtime(2)=wftime(indj+1) … … 175 197 if (lwindinterv.gt.idiffmax) nstop=3 176 198 199 call verboutput(verbosity,1,' leaving getfields') 200 177 201 end subroutine getfields -
branches/petra/src/gridcheck.f90
r27 r37 32 32 ! LAST UPDATE: 1997-10-10 * 33 33 ! * 34 !********************************************************************** 35 ! * 34 36 ! Update: 1999-02-08, global fields allowed, A. Stohl* 35 37 ! CHANGE: 11/01/2008, Harald Sodemann, GRIB1/2 input with * … … 37 39 ! CHANGE: 03/12/2008, Harald Sodemann, update to f90 with * 38 40 ! ECMWF grib_api * 39 ! * 41 ! 42 ! 3/2015, PS: add GRIB2 coding for deta/dt dp/deta used by ZAMG 43 ! 44 ! 40 45 !********************************************************************** 41 46 ! * … … 144 149 if (gribVer.eq.1) then ! GRIB Edition 1 145 150 146 !print*,'GRiB Edition 1' 147 !read the grib2 identifiers 148 call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret) 149 call grib_check(iret,gribFunction,gribErrorMsg) 150 call grib_get_int(igrib,'level',isec1(8),iret) 151 call grib_check(iret,gribFunction,gribErrorMsg) 152 153 !change code for etadot to code for omega 154 if (isec1(6).eq.77) then 155 isec1(6)=135 156 endif 157 158 !print*,isec1(6),isec1(8) 159 160 else 161 162 !print*,'GRiB Edition 2' 163 !read the grib2 identifiers 164 call grib_get_int(igrib,'discipline',discipl,iret) 165 call grib_check(iret,gribFunction,gribErrorMsg) 166 call grib_get_int(igrib,'parameterCategory',parCat,iret) 167 call grib_check(iret,gribFunction,gribErrorMsg) 168 call grib_get_int(igrib,'parameterNumber',parNum,iret) 169 call grib_check(iret,gribFunction,gribErrorMsg) 170 call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret) 171 call grib_check(iret,gribFunction,gribErrorMsg) 172 call grib_get_int(igrib,'level',valSurf,iret) 173 call grib_check(iret,gribFunction,gribErrorMsg) 174 call grib_get_int(igrib,'paramId',parId,iret) 175 call grib_check(iret,gribFunction,gribErrorMsg) 176 177 !print*,discipl,parCat,parNum,typSurf,valSurf 178 179 !convert to grib1 identifiers 180 isec1(6)=-1 181 isec1(7)=-1 182 isec1(8)=-1 183 isec1(8)=valSurf ! level 184 if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T 185 isec1(6)=130 ! indicatorOfParameter 186 elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.105)) then ! U 187 isec1(6)=131 ! indicatorOfParameter 188 elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.105)) then ! V 189 isec1(6)=132 ! indicatorOfParameter 190 elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q 191 isec1(6)=133 ! indicatorOfParameter 192 elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP 193 isec1(6)=134 ! indicatorOfParameter 194 elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot 195 isec1(6)=135 ! indicatorOfParameter 196 elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually eta dot 197 isec1(6)=135 ! indicatorOfParameter 198 elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP 199 isec1(6)=151 ! indicatorOfParameter 200 elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U 201 isec1(6)=165 ! indicatorOfParameter 202 elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V 203 isec1(6)=166 ! indicatorOfParameter 204 elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T 205 isec1(6)=167 ! indicatorOfParameter 206 elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D 207 isec1(6)=168 ! indicatorOfParameter 208 elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD 209 isec1(6)=141 ! indicatorOfParameter 210 elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC 211 isec1(6)=164 ! indicatorOfParameter 212 elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP 213 isec1(6)=142 ! indicatorOfParameter 214 elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP 215 isec1(6)=143 ! indicatorOfParameter 216 elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF 217 isec1(6)=146 ! indicatorOfParameter 218 elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR 219 isec1(6)=176 ! indicatorOfParameter 220 elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS 221 isec1(6)=180 ! indicatorOfParameter 222 elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS 223 isec1(6)=181 ! indicatorOfParameter 224 elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO 225 isec1(6)=129 ! indicatorOfParameter 226 elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO 227 isec1(6)=160 ! indicatorOfParameter 228 elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. & 229 (typSurf.eq.1)) then ! LSM 230 isec1(6)=172 ! indicatorOfParameter 231 else 232 print*,'***ERROR: undefined GRiB2 message found!',discipl, & 233 parCat,parNum,typSurf 234 endif 235 if(parId .ne. isec1(6) .and. parId .ne. 77) then 236 write(*,*) 'parId',parId, 'isec1(6)',isec1(6) 237 ! stop 238 endif 239 240 endif 151 !print*,'GRiB Edition 1' 152 !read the grib2 identifiers 153 call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret) 154 call grib_check(iret,gribFunction,gribErrorMsg) 155 call grib_get_int(igrib,'level',isec1(8),iret) 156 call grib_check(iret,gribFunction,gribErrorMsg) 157 158 !change code for etadot to code for omega 159 if (isec1(6).eq.77) then 160 isec1(6)=135 161 endif 162 163 !print*,isec1(6),isec1(8) 164 165 else ! GRIB 2 166 167 !print*,'GRiB Edition 2' 168 ! read the grib2 identifiers 169 170 call grib_get_int(igrib,'discipline',discipl,iret) 171 call grib_check(iret,gribFunction,gribErrorMsg) 172 call grib_get_int(igrib,'parameterCategory',parCat,iret) 173 call grib_check(iret,gribFunction,gribErrorMsg) 174 call grib_get_int(igrib,'parameterNumber',parNum,iret) 175 call grib_check(iret,gribFunction,gribErrorMsg) 176 call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret) 177 call grib_check(iret,gribFunction,gribErrorMsg) 178 call grib_get_int(igrib,'level',valSurf,iret) 179 call grib_check(iret,gribFunction,gribErrorMsg) 180 call grib_get_int(igrib,'paramId',parId,iret) 181 call grib_check(iret,gribFunction,gribErrorMsg) 182 183 !print*,discipl,parCat,parNum,typSurf,valSurf 184 185 ! convert to grib1 identifiers 186 187 isec1(6)=-1 188 isec1(7)=-1 189 isec1(8)=-1 190 isec1(8)=valSurf ! level 191 if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T 192 isec1(6)=130 ! indicatorOfParameter 193 elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.105)) then ! U 194 isec1(6)=131 ! indicatorOfParameter 195 elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.105)) then ! V 196 isec1(6)=132 ! indicatorOfParameter 197 elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q 198 isec1(6)=133 ! indicatorOfParameter 199 elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP 200 isec1(6)=134 ! indicatorOfParameter 201 elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually detadtdpdeta 202 isec1(6)=135 ! indicatorOfParameter 203 elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually detadtdpdeta 204 isec1(6)=135 ! indicatorOfParameter 205 elseif ((parCat.eq.2).and.(parNum.eq.8)) then ! Omega, actually detadtdpdeta 206 isec1(6)=135 ! indicatorOfParameter 207 elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP 208 isec1(6)=151 ! indicatorOfParameter 209 elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U 210 isec1(6)=165 ! indicatorOfParameter 211 elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V 212 isec1(6)=166 ! indicatorOfParameter 213 elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T 214 isec1(6)=167 ! indicatorOfParameter 215 elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D 216 isec1(6)=168 ! indicatorOfParameter 217 elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD 218 isec1(6)=141 ! indicatorOfParameter 219 elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC 220 isec1(6)=164 ! indicatorOfParameter 221 elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP 222 isec1(6)=142 ! indicatorOfParameter 223 elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP 224 isec1(6)=143 ! indicatorOfParameter 225 elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF 226 isec1(6)=146 ! indicatorOfParameter 227 elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR 228 isec1(6)=176 ! indicatorOfParameter 229 elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS 230 isec1(6)=180 ! indicatorOfParameter 231 elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS 232 isec1(6)=181 ! indicatorOfParameter 233 elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO 234 isec1(6)=129 ! indicatorOfParameter 235 elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO 236 isec1(6)=160 ! indicatorOfParameter 237 elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. & 238 (typSurf.eq.1)) then ! LSM 239 isec1(6)=172 ! indicatorOfParameter 240 else 241 print*,'***ERROR: undefined GRiB2 message found!',discipl, & 242 parCat,parNum,typSurf 243 endif 244 if(parId .ne. isec1(6) .and. parId .ne. 77) then 245 write(*,*) 'parId',parId, 'isec1(6)',isec1(6) 246 ! stop 247 endif 248 249 endif ! GRIB 1 / GRIB 2 cases 241 250 242 251 !get the size and data of the values array -
branches/petra/src/interpol_rain.f90
r24 r37 1 1 !********************************************************************** 2 ! Copyright 1998 ,1999,2000,2001,2002,2005,2007,2008,2009,2010*2 ! Copyright 1998-2015 * 3 3 ! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, * 4 4 ! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann * … … 20 20 !********************************************************************** 21 21 22 subroutine interpol_rain(yy1,yy2,yy3,nxmax,nymax,nzmax,nx, & 23 ny,memind,xt,yt,level,itime1,itime2,itime,yint1,yint2,yint3) 22 subroutine interpol_rain(yy1,yy2,yy3,nxmax,nymax,nzmax,nx,ny, & 23 memind,xt,yt,level,itime1,itime2,itime,yint1,yint2,yint3, & 24 icbot,ictop) 24 25 ! i i i i i i i 25 26 !i i i i i i i i o o o … … 40 41 ! * 41 42 ! 30 August 1996 * 43 !**************************************************************************** 44 ! Petra Seibert, 2011/2012-2015: 45 ! Add interpolation of cloud bottom and cloud thickness 46 ! for fix to SE's new wet scavenging scheme 42 47 ! * 43 48 !**************************************************************************** … … 66 71 ! * 67 72 !**************************************************************************** 73 74 use com_mod, only: icloudbot, icloudthck 75 use par_mod, only: icmv 68 76 69 77 implicit none 70 78 71 79 integer :: nx,ny,nxmax,nymax,nzmax,memind(2),m,ix,jy,ixp,jyp 72 integer :: itime,itime1,itime2,level,indexh 80 integer :: itime,itime1,itime2,level,indexh,ip1,ip2,ip3,ip4 81 integer :: icbot,icthck,ictop,ipsum 73 82 real :: yy1(0:nxmax-1,0:nymax-1,nzmax,2) 74 83 real :: yy2(0:nxmax-1,0:nymax-1,nzmax,2) 75 84 real :: yy3(0:nxmax-1,0:nymax-1,nzmax,2) 76 real :: ddx,ddy,rddx,rddy,dt1,dt2,dt,y1(2),y2(2),y3(2) 85 real :: ddx,ddy,rddx,rddy,dt1,dt2,dt,y1(2),y2(2),y3(2),cbot(2),cthck(2) 77 86 real :: xt,yt,yint1,yint2,yint3,p1,p2,p3,p4 78 79 87 80 88 … … 112 120 !*********************** 113 121 122 memind_loop: & 114 123 do m=1,2 124 115 125 indexh=memind(m) 116 126 … … 127 137 + p3*yy3(ix ,jyp,level,indexh) & 128 138 + p4*yy3(ixp,jyp,level,indexh) 129 end do 139 140 141 ! PS clouds: 142 ip1=1 143 ip2=1 144 ip3=1 145 ip4=1 146 if (icloudbot(ix ,jy ,indexh) .eq. icmv) ip1=0 147 if (icloudbot(ixp,jy ,indexh) .eq. icmv) ip2=0 148 if (icloudbot(ix ,jyp,indexh) .eq. icmv) ip3=0 149 if (icloudbot(ixp,jyp,indexh) .eq. icmv) ip4=0 150 ipsum = ip1+ip2+ip3+ip4 151 if (ipsum .eq. 0) then 152 cbot(m)=icmv 153 else 154 cbot(m)=(ip1*p1*icloudbot(ix ,jy ,indexh) & 155 + ip2*p2*icloudbot(ixp,jy ,indexh) & 156 + ip3*p3*icloudbot(ix ,jyp,indexh) & 157 + ip4*p4*icloudbot(ixp,jyp,indexh)) / ipsum 158 endif 159 160 ip1=1 161 ip2=1 162 ip3=1 163 ip4=1 164 if (icloudthck(ix ,jy ,indexh) .eq. icmv) ip1=0 165 if (icloudthck(ixp,jy ,indexh) .eq. icmv) ip2=0 166 if (icloudthck(ix ,jyp,indexh) .eq. icmv) ip3=0 167 if (icloudthck(ixp,jyp,indexh) .eq. icmv) ip4=0 168 ipsum = ip1+ip2+ip3+ip4 169 if (ipsum .eq. 0) then 170 cthck(m)=icmv 171 else 172 cthck(m)=(ip1*p1*icloudthck(ix ,jy ,indexh) & 173 + ip2*p2*icloudthck(ixp,jy ,indexh) & 174 + ip3*p3*icloudthck(ix ,jyp,indexh) & 175 + ip4*p4*icloudthck(ixp,jyp,indexh)) / ipsum 176 endif 177 ! PS end clouds 178 179 enddo memind_loop 130 180 131 181 … … 142 192 yint3=(y3(1)*dt2+y3(2)*dt1)/dt 143 193 194 ! PS clouds: 195 icbot = nint( (cbot(1)*dt2 + cbot(2)*dt1)/dt ) 196 if (nint(cbot(1)) .eq. icmv) icbot=cbot(2) 197 if (nint(cbot(2)) .eq. icmv) icbot=cbot(1) 198 199 icthck = nint( (cthck(1)*dt2 + cthck(2)*dt1)/dt ) 200 if (nint(cthck(1)) .eq. icmv) icthck=cthck(2) 201 if (nint(cthck(2)) .eq. icmv) icthck=cthck(1) 202 203 if (icbot .ne. icmv .and. icthck .ne. icmv) then 204 ictop = icbot + icthck ! convert cloud thickness to cloud top 205 else 206 icbot=icmv 207 ictop=icmv 208 endif 209 ! PS end clouds 210 144 211 145 212 end subroutine interpol_rain -
branches/petra/src/interpol_rain_nests.f90
r4 r37 1 1 !********************************************************************** 2 ! Copyright 1998 ,1999,2000,2001,2002,2005,2007,2008,2009,2010*2 ! Copyright 1998-2015 * 3 3 ! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, * 4 4 ! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann * … … 21 21 22 22 subroutine interpol_rain_nests(yy1,yy2,yy3,nxmaxn,nymaxn,nzmax, & 23 maxnests,ngrid,nxn,nyn,memind,xt,yt,level,itime1,itime2,itime, & 24 yint1,yint2,yint3) 23 maxnests,ngrid,nxn,nyn, & 24 memind,xt,yt,level,itime1,itime2,itime,yint1,yint2,yint3, & 25 icbot,ictop) 25 26 ! i i i i i i 26 27 ! i i i i i i i i i i i … … 44 45 ! * 45 46 ! 15 March 2000 * 47 ! * 48 !**************************************************************************** 49 ! Petra Seibert, 2011/2012-2015: 50 ! Add interpolation of cloud bottom and cloud thickness 51 ! for fix to SE's new wet scavenging scheme 46 52 ! * 47 53 !**************************************************************************** … … 71 77 !**************************************************************************** 72 78 79 use com_mod, only: icloudbotn, icloudthckn 80 use par_mod, only: icmv 81 73 82 implicit none 74 83 75 84 integer :: maxnests,ngrid 76 integer :: nxn(maxnests),nyn(maxnests),nxmaxn,nymaxn,nzmax,memind(2) 77 integer :: m,ix,jy,ixp,jyp,itime,itime1,itime2,level,indexh 85 integer :: nxn(maxnests),nyn(maxnests),nxmaxn,nymaxn,nzmax,memind(2),m 86 integer :: ix,jy,ixp,jyp 87 integer :: itime,itime1,itime2,level,indexh,ip1,ip2,ip3,ip4 88 integer :: icbot,ictop,icthck,ipsum 78 89 real :: yy1(0:nxmaxn-1,0:nymaxn-1,nzmax,2,maxnests) 79 90 real :: yy2(0:nxmaxn-1,0:nymaxn-1,nzmax,2,maxnests) 80 91 real :: yy3(0:nxmaxn-1,0:nymaxn-1,nzmax,2,maxnests) 81 real :: ddx,ddy,rddx,rddy,dt1,dt2,dt,y1(2),y2(2),y3(2) 92 real :: ddx,ddy,rddx,rddy,dt1,dt2,dt,y1(2),y2(2),y3(2),cbot(2),cthck(2) 82 93 real :: xt,yt,yint1,yint2,yint3,p1,p2,p3,p4 83 84 85 94 86 95 ! If point at border of grid -> small displacement into grid 87 96 !*********************************************************** 88 97 89 if (xt.ge.(real(nxn(ngrid)-1)-0.0001)) & 90 xt=real(nxn(ngrid)-1)-0.0001 91 if (yt.ge.(real(nyn(ngrid)-1)-0.0001)) & 92 yt=real(nyn(ngrid)-1)-0.0001 98 if (xt.ge.(real(nxn(ngrid)-1)-0.0001)) xt=real(nxn(ngrid)-1)-0.0001 99 if (yt.ge.(real(nyn(ngrid)-1)-0.0001)) yt=real(nyn(ngrid)-1)-0.0001 93 100 94 101 … … 119 126 !*********************** 120 127 128 memind_loop: & 121 129 do m=1,2 130 122 131 indexh=memind(m) 123 132 … … 134 143 + p3*yy3(ix ,jyp,level,indexh,ngrid) & 135 144 + p4*yy3(ixp,jyp,level,indexh,ngrid) 136 end do 145 146 ! PS clouds: 147 ip1=1 148 ip2=1 149 ip3=1 150 ip4=1 151 if (icloudbotn(ix ,jy ,indexh,ngrid) .eq. icmv) ip1=0 152 if (icloudbotn(ixp,jy ,indexh,ngrid) .eq. icmv) ip2=0 153 if (icloudbotn(ix ,jyp,indexh,ngrid) .eq. icmv) ip3=0 154 if (icloudbotn(ixp,jyp,indexh,ngrid) .eq. icmv) ip4=0 155 ipsum = ip1+ip2+ip3+ip4 156 if (ipsum .eq. 0) then 157 cbot(m)=icmv 158 else 159 cbot(m)=(ip1*p1*icloudbotn(ix ,jy ,indexh,ngrid) & 160 + ip2*p2*icloudbotn(ixp,jy ,indexh,ngrid) & 161 + ip3*p3*icloudbotn(ix ,jyp,indexh,ngrid) & 162 + ip4*p4*icloudbotn(ixp,jyp,indexh,ngrid)) / ipsum 163 endif 164 165 ip1=1 166 ip2=1 167 ip3=1 168 ip4=1 169 if (icloudthckn(ix ,jy ,indexh,ngrid) .eq. icmv) ip1=0 170 if (icloudthckn(ixp,jy ,indexh,ngrid) .eq. icmv) ip2=0 171 if (icloudthckn(ix ,jyp,indexh,ngrid) .eq. icmv) ip3=0 172 if (icloudthckn(ixp,jyp,indexh,ngrid) .eq. icmv) ip4=0 173 ipsum = ip1+ip2+ip3+ip4 174 if (ipsum .eq. 0) then 175 cthck(m)=icmv 176 else 177 cthck(m)=(ip1*p1*icloudthckn(ix ,jy ,indexh,ngrid) & 178 + ip2*p2*icloudthckn(ixp,jy ,indexh,ngrid) & 179 + ip3*p3*icloudthckn(ix ,jyp,indexh,ngrid) & 180 + ip4*p4*icloudthckn(ixp,jyp,indexh,ngrid)) / ipsum 181 endif 182 ! PS end clouds 183 184 enddo memind_loop 137 185 138 186 … … 150 198 151 199 200 ! PS clouds: 201 icbot = nint( (cbot(1)*dt2 + cbot(2)*dt1)/dt ) 202 if (nint(cbot(1)) .eq. icmv) icbot=cbot(2) 203 if (nint(cbot(2)) .eq. icmv) icbot=cbot(1) 204 205 icthck = nint( (cthck(1)*dt2 + cthck(2)*dt1)/dt ) 206 if (nint(cthck(1)) .eq. icmv) icthck=cthck(2) 207 if (nint(cthck(2)) .eq. icmv) icthck=cthck(1) 208 209 if (icbot .ne. icmv .and. icthck .ne. icmv) then 210 ictop = icbot + icthck ! convert cloud thickness to cloud top 211 else 212 icbot=icmv 213 ictop=icmv 214 endif 215 ! PS end clouds 216 217 152 218 end subroutine interpol_rain_nests -
branches/petra/src/par_mod.f90
r36 r37 103 103 104 104 integer,parameter :: idiffnorm=10800, idiffmax=2*idiffnorm, minstep=1 105 integer,parameter :: itagekernmin=10800 105 106 106 107 ! idiffnorm [s] normal time interval between two wind fields 107 108 ! idiffmax [s] maximum time interval between two wind fields 108 109 ! minstep [s] minimum time step to be used within FLEXPART 109 110 ! itagekernmin [s] minimum particle age [s] for using output kernel 110 111 111 112 !***************************************************************** … … 265 266 266 267 !****************************************************** 267 ! integer code for missing values, used in wet scavenging (PS, 2012)268 ! for use in wet scavenging (PS, 2012, 2015) 268 269 !****************************************************** 269 270 270 ! integer icmv 271 ! parameter(icmv=-9999) 272 integer,parameter :: icmv=-9999 271 integer,parameter :: icmv=-9999 ! integer code for missing values 272 ! some cloud heights used when trying to construct reasonable cloud: 273 integer,parameter :: icloudtopconvmin = 6000 274 integer,parameter :: icloudtopmin = 3000 275 integer,parameter :: icloudbot1 = 500, icloudtop1 = 8000 276 integer,parameter :: icloudbot2 = 0, icloudtop2 = 10000 277 278 real, parameter :: precmin = 0.002 ! minimum prec in mm/h for cloud diagnostics 279 real, parameter :: rhmininit = 0.90 ! standard condition for presence of clouds 280 ! PS note that original by Sabine Eckhart was 80% 281 ! PS however, for T<-20 C we consider saturation over ice 282 ! PS so I think 90% should be enough 283 real, parameter :: rhminx = 0.30 ! extreme condition for presence of clouds 273 284 274 285 ! Parameters for testing -
branches/petra/src/readcommand.f90
r36 r37 235 235 read(unitcommand,*,iostat=icmdstat) linit_cond 236 236 if (icmdstat .gt. 0) & 237 print*, 'readcommand: linit_cond not read properly',icmdstat,linit_cond 238 if (old) call skplin(3,unitcommand) 239 read(unitcommand,*,iostat=icmdstat) surf_only 237 print*, 'readcommand: linit_cond not present or no read properly,& 238 & - check values: ',icmdstat,linit_cond 239 if (old) call skplin(3,unitcommand) 240 read(unitcommand,*,iostat=icmdstat) icmdstat,surf_only 240 241 if (icmdstat .gt. 0) & 241 print*, 'readcommand: linit_cond not read properly',icmdstat,surf_only 242 if (old) call skplin(3,unitcommand) 243 read(unitcommand,*,iostat=icmdstat) ldep_incr 242 print*, 'readcommand: surf_only not present or not read properly& 243 & - check values: ',icmdstat,surf_only 244 if (old) call skplin(3,unitcommand) 245 read(unitcommand,*,iostat=icmdstat) icmdstat,ldep_incr 244 246 if (icmdstat .gt. 0) & 245 print*, 'readcommand: linit_cond not read properly',icmdstat, ldep_incr 247 print*, 'readcommand: ldep_incr not present or not read properly& 248 & - check values: ',icmdstat, ldep_incr 246 249 close(unitcommand) 247 250 … … 416 419 417 420 418 ! For backward runs one releasefield for all releases makes no sense, 419 ! For quasilag and domainfill ioutputforechrelease is forbidden 420 !***************************************************************************** 421 422 if ((ldirect .lt. 0) .and. (ioutputforeachrelease .eq. 0)) then 423 write(*,*) '#### FLEXPART MODEL ERROR! FILE COMMAND: ####' 424 write(*,*) '#### FOR BACKWARD RUNS, IOUTPUTFOREACHRLEASE ####' 425 write(*,*) '#### MUST BE SET TO ONE! ####' 426 stop 427 endif 428 429 430 ! For backward runs one releasefield for all releases makes no sense, 421 ! For domainfill runs one releasefield for all releases makes no sense, 431 422 ! and is "forbidden" 432 423 !***************************************************************************** -
branches/petra/src/readspecies.f90
r36 r37 20 20 !********************************************************************** 21 21 22 subroutine readspecies(id_spec, pos_spec)22 subroutine readspecies(id_spec,id_pos) 23 23 24 24 !***************************************************************************** 25 25 ! * 26 26 ! This routine reads names and physical constants of chemical species/ * 27 ! radionuclides given in the parameter pos_spec*27 ! radionuclides given in the parameter id_pos * 28 28 ! * 29 29 ! Author: A. Stohl * … … 32 32 ! 33 33 ! Changes: * 34 ! * 34 35 ! N. Kristiansen, 31.01.2013: Including parameters for in-cloud scavenging * 35 36 ! * 36 37 ! HSO, 13 August 2013: added optional namelist input 37 38 ! PS, 2/2015: access= -> position= 39 ! PS, 4/2015: quick fix for bug associated with end=22 38 40 ! * 39 41 !***************************************************************************** … … 60 62 implicit none 61 63 62 integer :: i, pos_spec,j64 integer :: i, id_pos,j 63 65 integer :: idow,ihour,id_spec 64 66 character(len=3) :: aspecnumb … … 69 71 real :: pdsigma, pdryvel, pweightmolar, pohreact, pspec_ass, pkao 70 72 real :: pweta_in, pwetb_in, pwetc_in, pwetd_in 71 integer :: readerror73 integer :: ios 72 74 73 75 ! declare namelist … … 100 102 ! Open the SPECIES file and read species names and properties 101 103 !************************************************************ 102 specnum( pos_spec)=id_spec103 write(aspecnumb,'(i3.3)') specnum( pos_spec)104 specnum(id_pos)=id_spec 105 write(aspecnumb,'(i3.3)') specnum(id_pos) 104 106 open(unitspecies,file=path(1)(1:length(1))//'SPECIES/SPECIES_'//aspecnumb,status='old',form='formatted',err=998) 105 !write(*,*) 'reading SPECIES',specnum( pos_spec)107 !write(*,*) 'reading SPECIES',specnum(id_pos) 106 108 107 109 ASSSPEC=.FALSE. 108 110 109 111 ! try namelist input 110 read(unitspecies,species_params,iostat= readerror)112 read(unitspecies,species_params,iostat=ios) 111 113 close(unitspecies) 112 114 113 if ((pweightmolar.eq.-789.0).or.(readerror.ne.0)) then ! no namelist found 114 115 readerror=1 116 117 open(unitspecies,file=path(1)(1:length(1))//'SPECIES/SPECIES_'//aspecnumb,status='old',err=998) 115 if ((pweightmolar.eq.-789.0).or.(ios.ne.0)) then ! no namelist found 116 117 open(unitspecies,file=path(1) (1:length(1))//'SPECIES/SPECIES_'//aspecnumb, status='old', err=998) 118 118 119 119 do i=1,6 120 120 read(unitspecies,*) 121 121 end do 122 123 read(unitspecies,'(a10)',end=22) species(pos_spec) 124 ! write(*,*) species(pos_spec) 125 read(unitspecies,'(f18.1)',end=22) decay(pos_spec) 126 ! write(*,*) decay(pos_spec) 127 read(unitspecies,'(e18.1)',end=22) weta(pos_spec) 128 ! write(*,*) weta(pos_spec) 129 read(unitspecies,'(f18.2)',end=22) wetb(pos_spec) 130 ! write(*,*) wetb(pos_spec) 122 read(unitspecies,'(a10)',end=20) species(id_pos) 123 ! write(*,*) species(id_pos) 124 read(unitspecies,'(f18.1)',end=20) decay(id_pos) 125 ! write(*,*) decay(id_pos) 126 read(unitspecies,'(e18.1)',end=20) weta(id_pos) 127 ! write(*,*) weta(id_pos) 128 read(unitspecies,'(f18.2)',end=20) wetb(id_pos) 129 ! write(*,*) wetb(id_pos) 131 130 132 131 !*** NIK 31.01.2013: including in-cloud scavening parameters 133 read(unitspecies,'(e18.1)',end=2 2) weta_in(pos_spec)134 ! write(*,*) weta_in( pos_spec)135 read(unitspecies,'(f18.2)',end=2 2) wetb_in(pos_spec)136 ! write(*,*) wetb_in( pos_spec)137 read(unitspecies,'(f18.2)',end=2 2) wetc_in(pos_spec)138 ! write(*,*) wetc_in( pos_spec)139 read(unitspecies,'(f18.2)',end=2 2) wetd_in(pos_spec)140 ! write(*,*) wetd_in( pos_spec)141 142 read(unitspecies,'(f18.1)',end=2 2) reldiff(pos_spec)143 ! write(*,*) reldiff( pos_spec)144 read(unitspecies,'(e18.1)',end=2 2) henry(pos_spec)145 ! write(*,*) henry( pos_spec)146 read(unitspecies,'(f18.1)',end=2 2) f0(pos_spec)147 ! write(*,*) f0( pos_spec)148 read(unitspecies,'(e18.1)',end=2 2) density(pos_spec)149 ! write(*,*) density( pos_spec)150 read(unitspecies,'(e18.1)',end=2 2) dquer(pos_spec)151 ! write(*,*) dquer( pos_spec)152 read(unitspecies,'(e18.1)',end=2 2) dsigma(pos_spec)153 ! write(*,*) dsigma( pos_spec)154 read(unitspecies,'(f18.2)',end=2 2) dryvel(pos_spec)155 ! write(*,*) dryvel( pos_spec)156 read(unitspecies,'(f18.2)',end=2 2) weightmolar(pos_spec)157 ! write(*,*) weightmolar( pos_spec)158 read(unitspecies,'(e18.1)',end=2 2) ohreact(pos_spec)159 ! write(*,*) ohreact( pos_spec)160 read(unitspecies,'(i18)',end=2 2) spec_ass(pos_spec)161 ! write(*,*) spec_ass( pos_spec)162 read(unitspecies,'(f18.2)',end=2 2) kao(pos_spec)163 ! write(*,*) kao( pos_spec)164 165 pspecies=species( pos_spec)166 pdecay=decay( pos_spec)167 pweta=weta( pos_spec)168 pwetb=wetb( pos_spec)169 pweta_in=weta_in( pos_spec)170 pwetb_in=wetb_in( pos_spec)171 pwetc_in=wetc_in( pos_spec)172 pwetd_in=wetd_in( pos_spec)173 preldiff=reldiff( pos_spec)174 phenry=henry( pos_spec)175 pf0=f0( pos_spec)176 pdensity=density( pos_spec)177 pdquer=dquer( pos_spec)178 pdsigma=dsigma( pos_spec)179 pdryvel=dryvel( pos_spec)180 pweightmolar=weightmolar( pos_spec)181 pohreact=ohreact( pos_spec)182 pspec_ass=spec_ass( pos_spec)183 pkao=kao( pos_spec)132 read(unitspecies,'(e18.1)',end=20) weta_in(id_pos) 133 ! write(*,*) weta_in(id_pos) 134 read(unitspecies,'(f18.2)',end=20) wetb_in(id_pos) 135 ! write(*,*) wetb_in(id_pos) 136 read(unitspecies,'(f18.2)',end=20) wetc_in(id_pos) 137 ! write(*,*) wetc_in(id_pos) 138 read(unitspecies,'(f18.2)',end=20) wetd_in(id_pos) 139 ! write(*,*) wetd_in(id_pos) 140 141 read(unitspecies,'(f18.1)',end=20) reldiff(id_pos) 142 ! write(*,*) reldiff(id_pos) 143 read(unitspecies,'(e18.1)',end=20) henry(id_pos) 144 ! write(*,*) henry(id_pos) 145 read(unitspecies,'(f18.1)',end=20) f0(id_pos) 146 ! write(*,*) f0(id_pos) 147 read(unitspecies,'(e18.1)',end=20) density(id_pos) 148 ! write(*,*) density(id_pos) 149 read(unitspecies,'(e18.1)',end=20) dquer(id_pos) 150 ! write(*,*) dquer(id_pos) 151 read(unitspecies,'(e18.1)',end=20) dsigma(id_pos) 152 ! write(*,*) dsigma(id_pos) 153 read(unitspecies,'(f18.2)',end=20) dryvel(id_pos) 154 ! write(*,*) dryvel(id_pos) 155 read(unitspecies,'(f18.2)',end=20) weightmolar(id_pos) 156 ! write(*,*) weightmolar(id_pos) 157 read(unitspecies,'(e18.1)',end=20) ohreact(id_pos) 158 ! write(*,*) ohreact(id_pos) 159 read(unitspecies,'(i18)',end=20) spec_ass(id_pos) 160 ! write(*,*) spec_ass(id_pos) 161 read(unitspecies,'(f18.2)',end=20) kao(id_pos) 162 ! write(*,*) kao(id_pos) 163 164 pspecies=species(id_pos) 165 pdecay=decay(id_pos) 166 pweta=weta(id_pos) 167 pwetb=wetb(id_pos) 168 pweta_in=weta_in(id_pos) 169 pwetb_in=wetb_in(id_pos) 170 pwetc_in=wetc_in(id_pos) 171 pwetd_in=wetd_in(id_pos) 172 preldiff=reldiff(id_pos) 173 phenry=henry(id_pos) 174 pf0=f0(id_pos) 175 pdensity=density(id_pos) 176 pdquer=dquer(id_pos) 177 pdsigma=dsigma(id_pos) 178 pdryvel=dryvel(id_pos) 179 pweightmolar=weightmolar(id_pos) 180 pohreact=ohreact(id_pos) 181 pspec_ass=spec_ass(id_pos) 182 pkao=kao(id_pos) 184 183 185 184 else 186 185 187 species(pos_spec)=pspecies 188 decay(pos_spec)=pdecay 189 weta(pos_spec)=pweta 190 wetb(pos_spec)=pwetb 191 weta_in(pos_spec)=pweta_in 192 wetb_in(pos_spec)=pwetb_in 193 wetc_in(pos_spec)=pwetc_in 194 wetd_in(pos_spec)=pwetd_in 195 reldiff(pos_spec)=preldiff 196 henry(pos_spec)=phenry 197 f0(pos_spec)=pf0 198 density(pos_spec)=pdensity 199 dquer(pos_spec)=pdquer 200 dsigma(pos_spec)=pdsigma 201 dryvel(pos_spec)=pdryvel 202 weightmolar(pos_spec)=pweightmolar 203 ohreact(pos_spec)=pohreact 204 spec_ass(pos_spec)=pspec_ass 205 kao(pos_spec)=pkao 206 207 endif 208 209 i=pos_spec 210 211 if ((weta(pos_spec).gt.0).and.(henry(pos_spec).le.0)) then 212 if (dquer(pos_spec).le.0) goto 996 ! no particle, no henry set 213 endif 214 215 if (spec_ass(pos_spec).gt.0) then 186 species(id_pos)=pspecies 187 decay(id_pos)=pdecay 188 weta(id_pos)=pweta 189 wetb(id_pos)=pwetb 190 weta_in(id_pos)=pweta_in 191 wetb_in(id_pos)=pwetb_in 192 wetc_in(id_pos)=pwetc_in 193 wetd_in(id_pos)=pwetd_in 194 reldiff(id_pos)=preldiff 195 henry(id_pos)=phenry 196 f0(id_pos)=pf0 197 density(id_pos)=pdensity 198 dquer(id_pos)=pdquer 199 dsigma(id_pos)=pdsigma 200 dryvel(id_pos)=pdryvel 201 weightmolar(id_pos)=pweightmolar 202 ohreact(id_pos)=pohreact 203 spec_ass(id_pos)=pspec_ass 204 kao(id_pos)=pkao 205 206 endif 207 208 if ((weta(id_pos).gt.0).and.(henry(id_pos).le.0)) then 209 if (dquer(id_pos).le.0) goto 996 ! no particle, no henry set 210 endif 211 212 if (spec_ass(id_pos).gt.0) then 216 213 spec_found=.FALSE. 217 do j=1, pos_spec-1218 if (spec_ass( pos_spec).eq.specnum(j)) then219 spec_ass( pos_spec)=j214 do j=1,id_pos-1 215 if (spec_ass(id_pos).eq.specnum(j)) then 216 spec_ass(id_pos)=j 220 217 spec_found=.TRUE. 221 218 ASSSPEC=.TRUE. … … 227 224 endif 228 225 229 if (dsigma(i ).eq.1.) dsigma(i)=1.0001 ! avoid floating exception230 if (dsigma(i ).eq.0.) dsigma(i)=1.0001 ! avoid floating exception231 232 if ((reldiff(i ).gt.0.).and.(density(i).gt.0.)) then226 if (dsigma(id_pos).eq.1.) dsigma(id_pos)=1.0001 ! avoid floating exception 227 if (dsigma(id_pos).eq.0.) dsigma(id_pos)=1.0001 ! avoid floating exception 228 229 if ((reldiff(id_pos).gt.0.).and.(density(id_pos).gt.0.)) then 233 230 write(*,*) '#### FLEXPART MODEL ERROR! FILE "SPECIES" ####' 234 231 write(*,*) '#### IS CORRUPT. SPECIES CANNOT BE BOTH ####' … … 246 243 247 244 do j=1,24 ! initialize everything to no variation 248 area_hour(i ,j)=1.249 point_hour(i ,j)=1.245 area_hour(id_pos,j)=1. 246 point_hour(id_pos,j)=1. 250 247 end do 251 248 do j=1,7 252 area_dow(i ,j)=1.253 point_dow(i ,j)=1.249 area_dow(id_pos,j)=1. 250 point_dow(id_pos,j)=1. 254 251 end do 255 252 256 if (readerror.ne.0) then ! text format input 257 258 read(unitspecies,*,end=22) 253 if (ios.ne.0) then ! text format input 254 255 read(unitspecies,*,iostat=ios) 256 if (ios .ne. 0) goto 22 ! ifort does not like err= 259 257 do j=1,24 ! 24 hours, starting with 0-1 local time 260 read(unitspecies,*) ihour,area_hour(i ,j),point_hour(i,j)258 read(unitspecies,*) ihour,area_hour(id_pos,j),point_hour(id_pos,j) 261 259 end do 262 260 read(unitspecies,*) 263 261 do j=1,7 ! 7 days of the week, starting with Monday 264 read(unitspecies,*) idow,area_dow(i ,j),point_dow(i,j)262 read(unitspecies,*) idow,area_dow(id_pos,j),point_dow(id_pos,j) 265 263 end do 266 264 -
branches/petra/src/readwind.f90
r36 r37 40 40 ! Variables tth and qvh (on eta coordinates) in common block 41 41 ! 2/2015, PS: add missing paramter iret in call to grib subr 42 ! 3/2015, PS: add some verbosity messages 43 ! 3/2015, PS: add GRIB2 coding for deta/dt dp/deta used by ZAMG 42 44 ! 43 45 !********************************************************************** … … 113 115 ! OPENING OF DATA FILE (GRIB CODE) 114 116 ! 115 5 call grib_open_file(ifile,path(3)(1:length(3)) & 117 118 5 continue 119 call verboutput(verbosity,1,' call grib_open_file '//path(3)(1:length(3))& 120 //trim(wfname(indj))) 121 call grib_open_file(ifile,path(3)(1:length(3)) & 116 122 //trim(wfname(indj)),'r',iret) 117 123 if (iret.ne.GRIB_SUCCESS) then … … 127 133 ! GET NEXT FIELDS 128 134 ! 135 call verboutput(verbosity,2,' call grib_new_from_file') 129 136 call grib_new_from_file(ifile,igrib,iret) 130 137 if (iret.eq.GRIB_END_OF_FILE) then … … 135 142 136 143 !first see if we read GRIB1 or GRIB2 144 call verboutput(verbosity,2,' call grib_get_int') 137 145 call grib_get_int(igrib,'editionNumber',gribVer,iret) 146 call verboutput(verbosity,2,' call grib_check') 138 147 call grib_check(iret,gribFunction,gribErrorMsg) 139 148 140 149 if (gribVer.eq.1) then ! GRIB Edition 1 141 150 142 !print*,'GRiB Edition 1' 143 !read the grib2 identifiers 144 call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret) 145 call grib_check(iret,gribFunction,gribErrorMsg) 146 call grib_get_int(igrib,'level',isec1(8),iret) 147 call grib_check(iret,gribFunction,gribErrorMsg) 148 149 !change code for etadot to code for omega 150 if (isec1(6).eq.77) then 151 isec1(6)=135 152 endif 153 154 conversion_factor=1. 155 156 else 157 158 !print*,'GRiB Edition 2' 159 !read the grib2 identifiers 160 call grib_get_int(igrib,'discipline',discipl,iret) 161 call grib_check(iret,gribFunction,gribErrorMsg) 162 call grib_get_int(igrib,'parameterCategory',parCat,iret) 163 call grib_check(iret,gribFunction,gribErrorMsg) 164 call grib_get_int(igrib,'parameterNumber',parNum,iret) 165 call grib_check(iret,gribFunction,gribErrorMsg) 166 call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret) 167 call grib_check(iret,gribFunction,gribErrorMsg) 168 call grib_get_int(igrib,'level',valSurf,iret) 169 call grib_check(iret,gribFunction,gribErrorMsg) 170 call grib_get_int(igrib,'paramId',parId,iret) 171 call grib_check(iret,gribFunction,gribErrorMsg) 172 173 !print*,discipl,parCat,parNum,typSurf,valSurf 174 175 !convert to grib1 identifiers 176 isec1(6)=-1 177 isec1(7)=-1 178 isec1(8)=-1 179 isec1(8)=valSurf ! level 180 conversion_factor=1. 181 if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T 182 isec1(6)=130 ! indicatorOfParameter 183 elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.105)) then ! U 184 isec1(6)=131 ! indicatorOfParameter 185 elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.105)) then ! V 186 isec1(6)=132 ! indicatorOfParameter 187 elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q 188 isec1(6)=133 ! indicatorOfParameter 189 elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP 190 isec1(6)=134 ! indicatorOfParameter 191 elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot 192 isec1(6)=135 ! indicatorOfParameter 193 elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually eta dot 194 isec1(6)=135 ! indicatorOfParameter 195 elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP 196 isec1(6)=151 ! indicatorOfParameter 197 elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U 198 isec1(6)=165 ! indicatorOfParameter 199 elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V 200 isec1(6)=166 ! indicatorOfParameter 201 elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T 202 isec1(6)=167 ! indicatorOfParameter 203 elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D 204 isec1(6)=168 ! indicatorOfParameter 205 elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD 206 isec1(6)=141 ! indicatorOfParameter 207 conversion_factor=1000. 208 elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC 209 isec1(6)=164 ! indicatorOfParameter 210 elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP 211 isec1(6)=142 ! indicatorOfParameter 212 elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP 213 isec1(6)=143 ! indicatorOfParameter 214 conversion_factor=1000. 215 elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF 216 isec1(6)=146 ! indicatorOfParameter 217 elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR 218 isec1(6)=176 ! indicatorOfParameter 219 elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS 220 isec1(6)=180 ! indicatorOfParameter 221 elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS 222 isec1(6)=181 ! indicatorOfParameter 223 elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO 224 isec1(6)=129 ! indicatorOfParameter 225 elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO 226 isec1(6)=160 ! indicatorOfParameter 227 elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. & 228 (typSurf.eq.1)) then ! LSM 229 isec1(6)=172 ! indicatorOfParameter 230 else 231 print*,'***WARNING: undefined GRiB2 message found!',discipl, & 232 parCat,parNum,typSurf 233 endif 234 if(parId .ne. isec1(6) .and. parId .ne. 77) then 235 write(*,*) 'parId',parId, 'isec1(6)',isec1(6) 236 ! stop 237 endif 238 239 endif 151 !print*,'GRiB Edition 1' 152 153 ! read the grib2 identifiers 154 155 call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret) 156 call grib_check(iret,gribFunction,gribErrorMsg) 157 call grib_get_int(igrib,'level',isec1(8),iret) 158 call grib_check(iret,gribFunction,gribErrorMsg) 159 160 !change code for etadot to code for omega 161 if (isec1(6).eq.77) then 162 isec1(6)=135 163 endif 164 165 conversion_factor=1. 166 167 else ! GRIB 2 168 169 !print*,'GRiB Edition 2' 170 171 ! read the grib2 identifiers 172 173 call grib_get_int(igrib,'discipline',discipl,iret) 174 call grib_check(iret,gribFunction,gribErrorMsg) 175 call grib_get_int(igrib,'parameterCategory',parCat,iret) 176 call grib_check(iret,gribFunction,gribErrorMsg) 177 call grib_get_int(igrib,'parameterNumber',parNum,iret) 178 call grib_check(iret,gribFunction,gribErrorMsg) 179 call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret) 180 call grib_check(iret,gribFunction,gribErrorMsg) 181 call grib_get_int(igrib,'level',valSurf,iret) 182 call grib_check(iret,gribFunction,gribErrorMsg) 183 call grib_get_int(igrib,'paramId',parId,iret) 184 call grib_check(iret,gribFunction,gribErrorMsg) 185 186 !print*,discipl,parCat,parNum,typSurf,valSurf 187 188 ! convert to grib1 identifiers 189 190 isec1(6)=-1 191 isec1(7)=-1 192 isec1(8)=-1 193 isec1(8)=valSurf ! level 194 conversion_factor=1. 195 if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T 196 isec1(6)=130 ! indicatorOfParameter 197 elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.105)) then ! U 198 isec1(6)=131 ! indicatorOfParameter 199 elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.105)) then ! V 200 isec1(6)=132 ! indicatorOfParameter 201 elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q 202 isec1(6)=133 ! indicatorOfParameter 203 elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP 204 isec1(6)=134 ! indicatorOfParameter 205 elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually detadtdpdeta 206 isec1(6)=135 ! indicatorOfParameter 207 elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually detadtdpdeta 208 isec1(6)=135 ! indicatorOfParameter 209 elseif ((parCat.eq.2).and.(parNum.eq.8)) then ! Omega, actually detadtdpdeta 210 isec1(6)=135 ! indicatorOfParameter 211 elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP 212 isec1(6)=151 ! indicatorOfParameter 213 elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U 214 isec1(6)=165 ! indicatorOfParameter 215 elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V 216 isec1(6)=166 ! indicatorOfParameter 217 elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T 218 isec1(6)=167 ! indicatorOfParameter 219 elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D 220 isec1(6)=168 ! indicatorOfParameter 221 elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD 222 isec1(6)=141 ! indicatorOfParameter 223 conversion_factor=1000. 224 elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC 225 isec1(6)=164 ! indicatorOfParameter 226 elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP 227 isec1(6)=142 ! indicatorOfParameter 228 elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP 229 isec1(6)=143 ! indicatorOfParameter 230 conversion_factor=1000. 231 elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF 232 isec1(6)=146 ! indicatorOfParameter 233 elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR 234 isec1(6)=176 ! indicatorOfParameter 235 elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS 236 isec1(6)=180 ! indicatorOfParameter 237 elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS 238 isec1(6)=181 ! indicatorOfParameter 239 elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO 240 isec1(6)=129 ! indicatorOfParameter 241 elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO 242 isec1(6)=160 ! indicatorOfParameter 243 elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. & 244 (typSurf.eq.1)) then ! LSM 245 isec1(6)=172 ! indicatorOfParameter 246 else 247 print*,'***WARNING: undefined GRiB2 message found!',discipl, & 248 parCat,parNum,typSurf 249 endif 250 if(parId .ne. isec1(6) .and. parId .ne. 77) then 251 write(*,*) 'parId',parId, 'isec1(6)',isec1(6) 252 ! stop 253 endif 254 255 endif ! end GRIB 1 / GRIB 2 cases 240 256 241 257 !HSO get the size and data of the values array … … 351 367 end do 352 368 369 call verboutput(verbosity,2,' call grib_release') 353 370 call grib_release(igrib) 354 371 goto 10 !! READ NEXT LEVEL OR PARAMETER … … 380 397 381 398 if (xglobal) then 399 call verboutput(verbosity,2,' call shift_field_0') 382 400 call shift_field_0(ewss,nxfield,ny) 383 401 call shift_field_0(nsss,nxfield,ny) … … 429 447 ff10m= sqrt(u10(i,j,1,n)**2+v10(i,j,1,n)**2) 430 448 fflev1=sqrt(uuh(i,j,3)**2+vvh(i,j,3)**2) 449 call verboutput(verbosity,2,' call pbl_profile') 431 450 call pbl_profile(ps(i,j,1,n),td2(i,j,1,n),hlev1, & 432 451 tt2(i,j,1,n),tth(i,j,3,n),ff10m,fflev1, & … … 457 476 if(iwmax.ne.nwz) stop 'READWIND: NWZ NOT CONSISTENT' 458 477 478 call verboutput(verbosity,1,' leaving readwind') 459 479 return 480 460 481 888 write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD #### ' 461 482 write(*,*) ' #### ',wfname(indj),' #### ' -
branches/petra/src/readwind_nests.f90
r36 r37 34 34 ! * 35 35 !***************************************************************************** 36 ! Changes, Bernd C. Krueger, Feb. 2001: * 36 ! CHANGES: 37 ! 38 ! 2/2001, Bernc C. Krueger: * 37 39 ! Variables tthn and qvhn (on eta coordinates) in common block * 38 ! CHANGE: 11/01/2008, Harald Sodemann, GRIB1/2 input with ECMWF grib_api * 39 ! CHANGE: 03/12/2008, Harald Sodemann, update to f90 with ECMWF grib_api * 40 ! 2/2015, PS: add missing paramter iret in call to grib subr 40 ! 11/01/2008, Harald Sodemann, GRIB1/2 input with ECMWF grib_api * 41 ! 03/12/2008, Harald Sodemann, update to f90 with ECMWF grib_api * 42 ! 2/2015, PS: add missing paramter iret in call to grib subr 43 ! 3/2015, PS: add verbosity output 44 ! 3/2015, PS: add GRIB2 coding for deta/dt dp/deta used by ZAMG 41 45 ! 42 46 !***************************************************************************** … … 84 88 character(len=24) :: gribErrorMsg = 'Error reading grib file' 85 89 character(len=20) :: gribFunction = 'readwind_nests' 90 91 call verboutput(verbosity,1,' entered readwind_nests') 86 92 87 93 do l=1,numbnests … … 100 106 ! 101 107 102 5 call grib_open_file(ifile,path(numpath+2*(l-1)+1) & 108 5 continue 109 call verboutput(verbosity,1,' call grib_open_file '//& 110 path(numpath+2*(l-1)+1)(1:length(numpath+2*(l-1)+1))//trim(wfnamen(l,indj))) 111 call grib_open_file(ifile,path(numpath+2*(l-1)+1) & 103 112 (1:length(numpath+2*(l-1)+1))//trim(wfnamen(l,indj)),'r') 104 113 if (iret.ne.GRIB_SUCCESS) then … … 127 136 if (gribVer.eq.1) then ! GRIB Edition 1 128 137 129 !print*,'GRiB Edition 1' 130 !read the grib2 identifiers 131 call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret) 132 call grib_check(iret,gribFunction,gribErrorMsg) 133 call grib_get_int(igrib,'level',isec1(8),iret) 134 call grib_check(iret,gribFunction,gribErrorMsg) 135 136 !change code for etadot to code for omega 137 if (isec1(6).eq.77) then 138 isec1(6)=135 139 endif 140 141 conversion_factor=1. 142 138 !print*,'GRiB Edition 1' 139 140 ! read the grib2 identifiers 141 142 call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret) 143 call grib_check(iret,gribFunction,gribErrorMsg) 144 call grib_get_int(igrib,'level',isec1(8),iret) 145 call grib_check(iret,gribFunction,gribErrorMsg) 146 147 !change code for etadot to code for omega 148 if (isec1(6).eq.77) then 149 isec1(6)=135 150 endif 151 152 conversion_factor=1. 143 153 144 154 else 145 155 146 !print*,'GRiB Edition 2' 147 !read the grib2 identifiers 148 call grib_get_int(igrib,'discipline',discipl,iret) 149 call grib_check(iret,gribFunction,gribErrorMsg) 150 call grib_get_int(igrib,'parameterCategory',parCat,iret) 151 call grib_check(iret,gribFunction,gribErrorMsg) 152 call grib_get_int(igrib,'parameterNumber',parNum,iret) 153 call grib_check(iret,gribFunction,gribErrorMsg) 154 call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret) 155 call grib_check(iret,gribFunction,gribErrorMsg) 156 call grib_get_int(igrib,'level',valSurf,iret) 157 call grib_check(iret,gribFunction,gribErrorMsg) 158 call grib_get_int(igrib,'paramId',parId,iret) !added by mc to make it consisitent with new readwind.f90 159 call grib_check(iret,gribFunction,gribErrorMsg) !added by mc to make it consisitent with new readwind.f90 160 161 !print*,discipl,parCat,parNum,typSurf,valSurf 162 163 !convert to grib1 identifiers 164 isec1(6)=-1 165 isec1(7)=-1 166 isec1(8)=-1 167 isec1(8)=valSurf ! level 168 conversion_factor=1. 169 if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T 170 isec1(6)=130 ! indicatorOfParameter 171 elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.105)) then ! U 172 isec1(6)=131 ! indicatorOfParameter 173 elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.105)) then ! V 174 isec1(6)=132 ! indicatorOfParameter 175 elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q 176 isec1(6)=133 ! indicatorOfParameter 177 elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP 178 isec1(6)=134 ! indicatorOfParameter 179 elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot ! 180 isec1(6)=135 ! indicatorOfParameter 181 elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually eta dot !added by mc to make it consisitent with new readwind.f90 182 isec1(6)=135 ! indicatorOfParameter !added by mc to make it consisitent with new readwind.f90 183 elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP 184 isec1(6)=151 ! indicatorOfParameter 185 elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U 186 isec1(6)=165 ! indicatorOfParameter 187 elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V 188 isec1(6)=166 ! indicatorOfParameter 189 elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T 190 isec1(6)=167 ! indicatorOfParameter 191 elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D 192 isec1(6)=168 ! indicatorOfParameter 193 elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD 194 isec1(6)=141 ! indicatorOfParameter 195 conversion_factor=1000. !added by mc to make it consisitent with new readwind.f90 196 elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC !added by mc to make it consisitent with new readwind.f90 197 isec1(6)=164 ! indicatorOfParameter 198 elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP !added by mc to make it consisitent with new readwind.f90 199 isec1(6)=142 ! indicatorOfParameter 200 elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP 201 isec1(6)=143 ! indicatorOfParameter 202 conversion_factor=1000. !added by mc to make it consisitent with new readwind.f90 203 elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF 204 isec1(6)=146 ! indicatorOfParameter 205 elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR 206 isec1(6)=176 ! indicatorOfParameter 207 elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS !added by mc to make it consisitent with new readwind.f90 208 isec1(6)=180 ! indicatorOfParameter 209 elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS !added by mc to make it consisitent with new readwind.f90 210 isec1(6)=181 ! indicatorOfParameter 211 elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO 212 isec1(6)=129 ! indicatorOfParameter 213 elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO !added by mc to make it consisitent with new readwind.f90 214 isec1(6)=160 ! indicatorOfParameter 215 elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. & 216 (typSurf.eq.1)) then ! LSM 217 isec1(6)=172 ! indicatorOfParameter 218 else 219 print*,'***WARNING: undefined GRiB2 message found!',discipl, & 220 parCat,parNum,typSurf 221 endif 222 if(parId .ne. isec1(6) .and. parId .ne. 77) then !added by mc to make it consisitent with new readwind.f90 223 write(*,*) 'parId',parId, 'isec1(6)',isec1(6) ! 224 ! stop 225 endif 156 !print*,'GRiB Edition 2' 157 158 ! read the grib2 identifiers 159 160 call grib_get_int(igrib,'discipline',discipl,iret) 161 call grib_check(iret,gribFunction,gribErrorMsg) 162 call grib_get_int(igrib,'parameterCategory',parCat,iret) 163 call grib_check(iret,gribFunction,gribErrorMsg) 164 call grib_get_int(igrib,'parameterNumber',parNum,iret) 165 call grib_check(iret,gribFunction,gribErrorMsg) 166 call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret) 167 call grib_check(iret,gribFunction,gribErrorMsg) 168 call grib_get_int(igrib,'level',valSurf,iret) 169 call grib_check(iret,gribFunction,gribErrorMsg) 170 call grib_get_int(igrib,'paramId',parId,iret) !added by mc to make it consistent with new readwind.f90 171 call grib_check(iret,gribFunction,gribErrorMsg) !added by mc to make it consistent with new readwind.f90 172 173 !print*,discipl,parCat,parNum,typSurf,valSurf 174 175 !convert to grib1 identifiers 176 isec1(6)=-1 177 isec1(7)=-1 178 isec1(8)=-1 179 isec1(8)=valSurf ! level 180 conversion_factor=1. 181 if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T 182 isec1(6)=130 ! indicatorOfParameter 183 elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.105)) then ! U 184 isec1(6)=131 ! indicatorOfParameter 185 elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.105)) then ! V 186 isec1(6)=132 ! indicatorOfParameter 187 elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q 188 isec1(6)=133 ! indicatorOfParameter 189 elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP 190 isec1(6)=134 ! indicatorOfParameter 191 elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually detadtdpdeta 192 isec1(6)=135 ! indicatorOfParameter 193 elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually detadtdpdeta 194 isec1(6)=135 ! indicatorOfParameter !added by mc to make it consistent with new readwind.f90 195 elseif ((parCat.eq.2).and.(parNum.eq.8)) then ! Omega, actually detadtdpdeta 196 isec1(6)=135 ! indicatorOfParameter !added by mc to make it consistent with new readwind.f90 197 elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP 198 isec1(6)=151 ! indicatorOfParameter 199 elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U 200 isec1(6)=165 ! indicatorOfParameter 201 elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V 202 isec1(6)=166 ! indicatorOfParameter 203 elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T 204 isec1(6)=167 ! indicatorOfParameter 205 elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D 206 isec1(6)=168 ! indicatorOfParameter 207 elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD 208 isec1(6)=141 ! indicatorOfParameter 209 conversion_factor=1000. !added by mc to make it consistent with new readwind.f90 210 elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC !added by mc to make it consistent with new readwind.f90 211 isec1(6)=164 ! indicatorOfParameter 212 elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP !added by mc to make it consistent with new readwind.f90 213 isec1(6)=142 ! indicatorOfParameter 214 elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP 215 isec1(6)=143 ! indicatorOfParameter 216 conversion_factor=1000. !added by mc to make it consistent with new readwind.f90 217 elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF 218 isec1(6)=146 ! indicatorOfParameter 219 elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR 220 isec1(6)=176 ! indicatorOfParameter 221 elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS !added by mc to make it consistent with new readwind.f90 222 isec1(6)=180 ! indicatorOfParameter 223 elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS !added by mc to make it consistent with new readwind.f90 224 isec1(6)=181 ! indicatorOfParameter 225 elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO 226 isec1(6)=129 ! indicatorOfParameter 227 elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO !added by mc to make it consistent with new readwind.f90 228 isec1(6)=160 ! indicatorOfParameter 229 elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. & 230 (typSurf.eq.1)) then ! LSM 231 isec1(6)=172 ! indicatorOfParameter 232 else 233 print*,'***WARNING: undefined GRiB2 message found!',discipl, & 234 parCat,parNum,typSurf 235 endif 236 if(parId .ne. isec1(6) .and. parId .ne. 77) then !added by mc to make it consistent with new readwind.f90 237 write(*,*) 'parId',parId, 'isec1(6)',isec1(6) ! 238 ! stop 239 endif 226 240 227 241 endif … … 246 260 ! CHECK GRID SPECIFICATIONS 247 261 if(isec2(2).ne.nxn(l)) stop & 248 'READWIND: NX NOT CONSISTENT FOR A NESTING LEVEL'262 'READWIND: NX NOT CONSISTENT FOR A NESTING LEVEL' 249 263 if(isec2(3).ne.nyn(l)) stop & 250 'READWIND: NY NOT CONSISTENT FOR A NESTING LEVEL'251 if(isec2(12)/2-1.ne.nlev_ec) stop 'READWIND: VERTICAL DISCRET&252 IZATION NOT CONSISTENT FOR A NESTING LEVEL'264 'READWIND: NY NOT CONSISTENT FOR A NESTING LEVEL' 265 if(isec2(12)/2-1.ne.nlev_ec) stop & 266 'READWIND: VERTICAL DISCRETIZATION NOT CONSISTENT FOR A NESTING LEVEL' 253 267 endif ! ifield 254 268 255 269 !HSO get the second part of the grid dimensions only from GRiB1 messages 256 if (isec1(6) .eq. 167 .and. (gotGrid.eq.0)) then ! !added by mc to make it consis itent with new readwind.f90270 if (isec1(6) .eq. 167 .and. (gotGrid.eq.0)) then ! !added by mc to make it consistent with new readwind.f90 257 271 call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', & 258 272 xauxin,iret) … … 294 308 zsec4(nxn(l)*(nyn(l)-j-1)+i+1) 295 309 if(isec1(6).eq.141) sdn(i,j,1,n,l)= &!! SNOW DEPTH 296 zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/conversion_factor !added by mc to make it consis itent with new readwind.f90!310 zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/conversion_factor !added by mc to make it consistent with new readwind.f90! 297 311 if(isec1(6).eq.151) msln(i,j,1,n,l)= &!! SEA LEVEL PRESS. 298 312 zsec4(nxn(l)*(nyn(l)-j-1)+i+1) … … 312 326 endif 313 327 if(isec1(6).eq.143) then !! CONVECTIVE PREC. 314 convprecn(i,j,1,n,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/conversion_factor !added by mc to make it consis itent with new readwind.f90328 convprecn(i,j,1,n,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/conversion_factor !added by mc to make it consistent with new readwind.f90 315 329 if (convprecn(i,j,1,n,l).lt.0.) convprecn(i,j,1,n,l)=0. 316 330 endif … … 421 435 end do 422 436 437 call verboutput(verbosity,1,' leaving readwind_nests') 438 423 439 return 440 424 441 888 write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD #### ' 425 442 write(*,*) ' #### ',wfnamen(l,indj),' FOR NESTING LEVEL #### ' … … 431 448 write(*,*) ' #### ',wfnamen(l,indj),' #### ' 432 449 write(*,*) ' #### CANNOT BE OPENED FOR NESTING LEVEL ',l,'####' 450 call verboutput(verbosity,1,' leaving readwind') 433 451 434 452 end subroutine readwind_nests -
branches/petra/src/timemanager.f90
r30 r37 137 137 !write(*,45) itime,numpart,gridtotalunc,wetgridtotalunc,drygridtotalunc 138 138 write(*,46) float(itime)/3600,itime,numpart 139 if (verbosity.gt.0) then 140 write (*,*) 'timemanager> starting simulation' 141 if (verbosity.gt.1) then 142 CALL SYSTEM_CLOCK(count_clock) 143 WRITE(*,*) 'timemanager> SYSTEM CLOCK',(count_clock - count_clock0)/real(count_rate) 144 endif 139 call verboutput(verbosity,1,'timemanager> starting simulation') 140 if (verbosity.gt.1) then 141 CALL SYSTEM_CLOCK(count_clock) 142 WRITE(*,*) 'timemanager> SYSTEM CLOCK', (count_clock-count_clock0) / real(count_rate) 145 143 endif 146 144 -
branches/petra/src/verttransform.f90
r24 r37 1 1 !********************************************************************** 2 ! Copyright 1998 ,1999,2000,2001,2002,2005,2007,2008,2009,2010*2 ! Copyright 1998-2015 * 3 3 ! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, * 4 4 ! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann * … … 50 50 ! added the variable cloud for use with scavenging - descr. in com_mod 51 51 !***************************************************************************** 52 ! Petra Seibert, Feb 2015: Add quick fix from 2013 53 !***************************************************************************** 52 54 ! * 53 55 ! Variables: * 54 56 ! nx,ny,nz field dimensions in x,y and z direction * 55 ! clouds(0:nxmax,0:nymax,0:nzmax,2) cloud field for wet deposition *56 57 ! uu(0:nxmax,0:nymax,nzmax,2) wind components in x-direction [m/s] * 57 58 ! vv(0:nxmax,0:nymax,nzmax,2) wind components in y-direction [m/s] * … … 69 70 implicit none 70 71 71 integer :: ix,jy,kz,iz,n,kmin,kl,klp,ix1,jy1,ixp,jyp,ixm,jym 72 integer :: rain_cloud_above,kz_inv72 integer :: ix,jy,kz,iz,n,kmin,kl,klp,ix1,jy1,ixp,jyp,ixm,jym,kz_inv 73 integer :: icloudtop 73 74 real :: f_qvsat,pressure 74 real :: rh,lsp,convp 75 real :: rh,lsp,convp,prec,rhmin 75 76 real :: rhoh(nuvzmax),pinmconv(nzmax) 76 77 real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi … … 85 86 real,parameter :: const=r_air/ga 86 87 87 logical :: init = .true. 88 logical :: init = .true., lconvectprec, lsearch 88 89 89 90 … … 477 478 478 479 479 !write (*,*) 'initializing clouds, n:',n,nymin1,nxmin1,nz 480 ! create a cloud and rainout/washout field, clouds occur where rh>80% 481 ! total cloudheight is stored at level 0 480 ! write (*,*) 'diagnosing clouds, n:',n,nymin1,nxmin1,nz 481 jy_loop: & 482 482 do jy=0,nymin1 483 483 do ix=0,nxmin1 484 rain_cloud_above=0484 485 485 lsp=lsprec(ix,jy,1,n) 486 486 convp=convprec(ix,jy,1,n) 487 cloudsh(ix,jy,n)=0 488 do kz_inv=1,nz-1 489 kz=nz-kz_inv+1 490 pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n) 491 rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n)) 492 clouds(ix,jy,kz,n)=0 493 if (rh.gt.0.8) then ! in cloud 494 if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation 495 rain_cloud_above=1 496 cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ & 497 height(kz)-height(kz-1) 498 if (lsp.ge.convp) then 499 clouds(ix,jy,kz,n)=3 ! lsp dominated rainout 500 else 501 clouds(ix,jy,kz,n)=2 ! convp dominated rainout 502 endif 503 else ! no precipitation 504 clouds(ix,jy,kz,n)=1 ! cloud 487 prec=lsp+convp 488 if (lsp .gt. convp) then ! prectype='lsp' 489 lconvectprec = .false. 490 else ! prectype='cp ' 491 lconvectprec = .true. 492 endif 493 icloudbot(ix,jy,n)=icmv 494 icloudtop=icmv ! this is just a local variable 495 rhmin=rhmininit ! just initialise in a way that cond is true 496 lsearch=.true. 497 498 cloudsearch_loop: & 499 do while (rhmin .ge. rhminx .and. lsearch) 500 ! give up for < rhminx 501 502 kz_loop: & 503 do kz_inv=1,nz-1 504 kz=nz-kz_inv+1 505 pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n) 506 rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n)) 507 if (rh .gt. rhmin) then 508 if (icloudbot(ix,jy,n) .eq. icmv) then 509 icloudtop=nint(height(kz)) ! use int to save memory 505 510 endif 506 else ! no cloud 507 if (rain_cloud_above.eq.1) then ! scavenging 508 if (lsp.ge.convp) then 509 clouds(ix,jy,kz,n)=5 ! lsp dominated washout 510 else 511 clouds(ix,jy,kz,n)=4 ! convp dominated washout 512 endif 513 endif 514 endif 515 end do 516 end do 517 end do 511 icloudbot(ix,jy,n)=nint(height(kz)) 512 endif 513 end do kz_loop 514 515 ! PS try to get a cloud thicker than 50 m 516 ! PS if there is at least precmin mm/h 517 if (prec .gt. precmin .and. & 518 ( icloudbot(ix,jy,n) .eq. icmv .or. & 519 icloudtop-icloudbot(ix,jy,n) .lt. 50)) then 520 rhmin = rhmin - 0.05 521 else 522 lsearch = .false. 523 endif 524 525 enddo cloudsearch_loop 526 527 ! PS implement a rough fix for badly represented convection 528 ! PS is based on looking at a limited set of comparison data 529 if (lconvectprec .and. icloudtop .lt. icloudtopconvmin .or. & 530 icloudbot(ix,jy,n) .lt. icloudtopmin .and. prec .gt. precmin) then 531 if (convp .lt. 0.1) then 532 icloudbot(ix,jy,n) = icloudbot1 533 icloudtop = icloudtop1 534 else 535 icloudbot(ix,jy,n) = icloudbot2 536 icloudtop = icloudtop2 537 endif 538 endif 539 540 if (icloudtop .ne. icmv) then 541 icloudthck(ix,jy,n) = icloudtop-icloudbot(ix,jy,n) 542 else 543 icloudthck(ix,jy,n) = icmv ! no cloud found whatsoever 544 endif 545 546 ! PS get rid of too thin clouds 547 if (icloudthck(ix,jy,n) .lt. 50) then 548 icloudbot(ix,jy,n)=icmv 549 icloudthck(ix,jy,n)=icmv 550 endif 551 552 enddo 553 enddo jy_loop 518 554 519 555 end subroutine verttransform -
branches/petra/src/verttransform_nests.f90
r24 r37 1 1 !********************************************************************** 2 ! Copyright 1998 ,1999,2000,2001,2002,2005,2007,2008,2009,2010*2 ! Copyright 1998-2015 * 3 3 ! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, * 4 4 ! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann * … … 52 52 ! add the variable cloud for use with scavenging - descr. in com_mod 53 53 !***************************************************************************** 54 ! Petra Seibert, Feb 2015: Add quick fix from 2013 55 !***************************************************************************** 54 56 ! * 55 57 ! Variables: * … … 69 71 implicit none 70 72 71 integer :: ix,jy,kz,iz,n,l,kmin,kl,klp,ix1,jy1,ixp,jyp 72 integer :: rain_cloud_above,kz_inv 73 real :: f_qvsat,pressure,rh,lsp,convp 74 real :: wzlev(nwzmax),rhoh(nuvzmax),pinmconv(nzmax) 75 real :: uvwzlev(0:nxmaxn-1,0:nymaxn-1,nzmax) 76 real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi,cosf 73 integer :: ix,jy,kz,iz,n,l,kmin,kl,klp,ix1,jy1,ixp,jyp,kz_inv 74 integer :: icloudtop 75 real :: f_qvsat,pressure 76 real :: rh,lsp,convp,prec,rhmin 77 real :: rhoh(nuvzmax),pinmconv(nzmax) 78 real :: wzlev(nwzmax),uvwzlev(0:nxmaxn-1,0:nymaxn-1,nzmax) 79 real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi 77 80 real :: dzdx,dzdy 78 real :: dzdx1,dzdx2,dzdy1,dzdy2 81 real :: dzdx1,dzdx2,dzdy1,dzdy2,cosf 79 82 real :: uuhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests) 80 83 real :: vvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests) … … 83 86 real,parameter :: const=r_air/ga 84 87 88 logical :: lconvectprec, lsearch 85 89 86 90 ! Loop over all nests 87 91 !******************** 88 92 93 numbnest_loop: & 89 94 do l=1,numbnests 90 95 … … 143 148 144 149 145 ! Levels ,where u,v,t and q are given146 !*********************************** *150 ! Levels where u,v,t and q are given 151 !*********************************** 147 152 148 153 uun(ix,jy,1,n,l)=uuhn(ix,jy,1,l) … … 194 199 195 200 196 ! Levels ,where w is given197 !************************ *201 ! Levels where w is given 202 !************************ 198 203 199 204 wwn(ix,jy,1,n,l)=wwhn(ix,jy,1,l)*pinmconv(1) … … 275 280 276 281 end do 277 278 282 end do 279 283 end do 280 284 281 282 !write (*,*) 'initializing nested cloudsn, n:',n 283 ! create a cloud and rainout/washout field, cloudsn occur where rh>80% 285 ! write (*,*) 'diagnosing nested clouds, n:',n 286 jy_loop: & 284 287 do jy=0,nyn(l)-1 285 288 do ix=0,nxn(l)-1 286 rain_cloud_above=0 289 287 290 lsp=lsprecn(ix,jy,1,n,l) 288 291 convp=convprecn(ix,jy,1,n,l) 289 cloudsnh(ix,jy,n,l)=0 290 do kz_inv=1,nz-1 291 kz=nz-kz_inv+1 292 pressure=rhon(ix,jy,kz,n,l)*r_air*ttn(ix,jy,kz,n,l) 293 rh=qvn(ix,jy,kz,n,l)/f_qvsat(pressure,ttn(ix,jy,kz,n,l)) 294 cloudsn(ix,jy,kz,n,l)=0 295 if (rh.gt.0.8) then ! in cloud 296 if ((lsp.gt.0.01).or.(convp.gt.0.01)) then 297 rain_cloud_above=1 298 cloudsnh(ix,jy,n,l)=cloudsnh(ix,jy,n,l)+ & 299 height(kz)-height(kz-1) 300 if (lsp.ge.convp) then 301 cloudsn(ix,jy,kz,n,l)=3 ! lsp dominated rainout 302 else 303 cloudsn(ix,jy,kz,n,l)=2 ! convp dominated rainout 304 endif 305 else ! no precipitation 306 cloudsn(ix,jy,kz,n,l)=1 ! cloud 292 prec=lsp+convp 293 if (lsp .gt. convp) then ! prectype='lsp' 294 lconvectprec = .false. 295 else ! prectype='cp ' 296 lconvectprec = .true. 297 endif 298 icloudbotn(ix,jy,n,l)=icmv 299 icloudtop=icmv ! this is just a local variable 300 rhmin=rhmininit 301 lsearch=.true. 302 303 cloudsearch_loop: & 304 do while (rhmin .ge. rhminx .and. lsearch) 305 ! give up for < rhminx 306 307 kz_loop: & 308 do kz_inv=1,nz-1 309 kz=nz-kz_inv+1 310 pressure=rhon(ix,jy,kz,n,l)*r_air*ttn(ix,jy,kz,n,l) 311 rh=qvn(ix,jy,kz,n,l)/f_qvsat(pressure,ttn(ix,jy,kz,n,l)) 312 if (rh .gt. rhmin) then 313 if (icloudbotn(ix,jy,n,l) .eq. icmv) then 314 icloudtop=nint(height(kz)) ! use int to save memory 307 315 endif 308 else ! no cloud 309 if (rain_cloud_above.eq.1) then ! scavenging 310 if (lsp.ge.convp) then 311 cloudsn(ix,jy,kz,n,l)=5 ! lsp dominated washout 312 else 313 cloudsn(ix,jy,kz,n,l)=4 ! convp dominated washout 314 endif 315 endif 316 endif 317 end do 318 end do 319 end do 320 321 end do 316 icloudbotn(ix,jy,n,l)=nint(height(kz)) 317 endif 318 end do kz_loop 319 320 ! PS try to get a cloud thicker than 50 m 321 ! PS if there is at least precmin mm/h 322 if (prec .gt. precmin .and. & 323 ( icloudbotn(ix,jy,n,l) .eq. icmv .or. & 324 icloudtop-icloudbotn(ix,jy,n,l) .lt. 50)) then 325 rhmin = rhmin - 0.05 326 else 327 lsearch = .false. 328 endif 329 330 enddo cloudsearch_loop 331 332 ! PS implement a rough fix for badly represented convection 333 ! PS is based on looking at a limited set of comparison data 334 if (lconvectprec .and. icloudtop .lt. icloudtopconvmin .or. & 335 icloudbotn(ix,jy,n,l) .lt. icloudtopmin .and. prec .gt. precmin) then 336 if (convp .lt. 0.1) then 337 icloudbotn(ix,jy,n,l) = icloudbot1 338 icloudtop = icloudtop1 339 else 340 icloudbotn(ix,jy,n,l) = icloudbot2 341 icloudtop = icloudtop2 342 endif 343 endif 344 if (icloudtop .ne. icmv) then 345 icloudthckn(ix,jy,n,l) = icloudtop-icloudbotn(ix,jy,n,l) 346 else 347 icloudthckn(ix,jy,n,l) = icmv 348 endif 349 ! PS get rid of too thin clouds 350 if (icloudthck(ix,jy,n) .lt. 50) then 351 icloudbotn(ix,jy,n,l)=icmv 352 icloudthckn(ix,jy,n,l)=icmv 353 endif 354 355 enddo 356 enddo jy_loop 357 358 enddo numbnest_loop 322 359 323 360 end subroutine verttransform_nests -
branches/petra/src/wetdepo.f90
r24 r37 1 1 !********************************************************************** 2 ! Copyright 1998 ,1999,2000,2001,2002,2005,2007,2008,2009,2010*2 ! Copyright 1998-2015 * 3 3 ! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, * 4 4 ! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann * … … 39 39 ! Code may not be correct for decay of deposition! * 40 40 ! * 41 ! PS, 2/2015: implement wet depo quick fix from 2011/12 in this version 42 ! it implements interpolated and improved clouds 43 ! Also, certain deficiencies for thin clouds fixed also 44 ! Pass itage to wetdepokernel 45 ! * 41 46 !***************************************************************************** 42 47 ! * … … 71 76 72 77 integer :: jpart,itime,ltsample,loutnext,ldeltat,i,j,ix,jy 73 integer :: ngrid,itage,nage, hz,il,interp_time, n, clouds_v74 integer :: ks, kp 78 integer :: ngrid,itage,nage,kz,il,interp_time,n 79 integer :: ks, kp, n1,n2, icbot,ictop, indcloud 75 80 real :: S_i, act_temp, cl, cle ! in cloud scavenging 76 real :: clouds_h ! cloud height for the specific grid point 77 real :: xtn,ytn,lsp,convp,cc,grfraction,prec,wetscav 81 real :: xtn,ytn,lsp,convp,cc,grfraction,prec,wetscav,wetscavold,precsub,f 78 82 real :: wetdeposit(maxspec),restmass 79 83 real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled … … 97 101 !************************ 98 102 103 particle_loop: & 99 104 do jpart=1,numpart 100 105 101 106 if (itra1(jpart).eq.-999999999) goto 20 102 if (ldirect.eq.1)then107 if (ldirect.eq.1) then 103 108 if (itra1(jpart).gt.itime) goto 20 104 109 else 105 110 if (itra1(jpart).lt.itime) goto 20 106 111 endif 112 107 113 ! Determine age class of the particle 108 114 itage=abs(itra1(jpart)-itramem(jpart)) … … 123 129 goto 23 124 130 endif 125 end 131 enddo 126 132 23 continue 127 133 … … 131 137 132 138 if (ngrid.gt.0) then 133 xtn= (xtra1(jpart)-xln(ngrid))*xresoln(ngrid)134 ytn= (ytra1(jpart)-yln(ngrid))*yresoln(ngrid)139 xtn=real(xtra1(jpart)-xln(ngrid))*xresoln(ngrid) 140 ytn=real(ytra1(jpart)-yln(ngrid))*yresoln(ngrid) 135 141 ix=int(xtn) 136 142 jy=int(ytn) … … 148 154 149 155 if (ngrid.eq.0) then 150 call interpol_rain(lsprec,convprec,tcc,nxmax,nymax, 151 1,nx,ny,memind,real(xtra1(jpart)),real(ytra1(jpart)),1, &152 memtime(1),memtime(2),interp_time,lsp,convp,cc)156 call interpol_rain(lsprec,convprec,tcc,nxmax,nymax,1,nx,ny,memind,& 157 real(xtra1(jpart)),real(ytra1(jpart)), & 158 1,memtime(1),memtime(2),interp_time,lsp,convp,cc,icbot,ictop) 153 159 else 154 160 call interpol_rain_nests(lsprecn,convprecn,tccn, & 155 nxmaxn,nymaxn,1,maxnests,ngrid,nxn,nyn,memind,xtn,ytn,1, & 156 memtime(1),memtime(2),interp_time,lsp,convp,cc) 157 endif 158 159 if ((lsp.lt.0.01).and.(convp.lt.0.01)) goto 20 161 nxmaxn,nymaxn,1,maxnests,ngrid,nxn,nyn,memind,xtn,ytn, & 162 1,memtime(1),memtime(2),interp_time,lsp,convp,cc,icbot,ictop) 163 endif 164 165 ! PS 2012/2015: subtract a small value, eg 0.01 mm/h, 166 ! to remove spurious precip; replaces previous code 167 prec = lsp+convp 168 precsub = 0.01 169 if (prec .lt. precsub) then 170 goto 20 171 else 172 f = (prec-precsub)/prec 173 lsp = f*lsp 174 convp = f*convp 175 endif 176 160 177 ! get the level were the actual particle is in 161 178 do il=2,nz 162 179 if (height(il).gt.ztra1(jpart)) then 163 hz=il-1180 kz=il-1 164 181 goto 26 165 182 endif 166 end 183 enddo 167 184 26 continue 168 185 169 186 n=memind(2) 170 187 if (abs(memtime(1)-interp_time).lt.abs(memtime(2)-interp_time)) & 171 n=memind(1)188 n=memind(1) 172 189 173 190 ! if there is no precipitation or the particle is above the clouds no 174 191 ! scavenging is done 175 176 if (ngrid.eq.0) then 177 clouds_v=clouds(ix,jy,hz,n) 178 clouds_h=cloudsh(ix,jy,n) 179 else 180 clouds_v=cloudsn(ix,jy,hz,n,ngrid) 181 clouds_h=cloudsnh(ix,jy,n,ngrid) 182 endif 183 !write(*,*) 'there is 184 ! + precipitation',(clouds(ix,jy,ihz,n),ihz=1,20),lsp,convp,hz 185 if (clouds_v.le.1) goto 20 186 !write (*,*) 'there is scavenging' 192 ! PS: part of 2011/2012/2015 fix, replaces previous code 193 194 if (ztra1(jpart) .le. float(ictop)) then 195 if (ztra1(jpart) .gt. float(icbot)) then 196 indcloud = 2 ! in-cloud 197 else 198 indcloud = 1 ! below-cloud 199 endif 200 elseif (ictop .eq. icmv) then 201 indcloud = 0 ! no cloud found, use old scheme 202 else 203 goto 20 ! above cloud 204 endif 205 187 206 188 207 ! 1) Parameterization of the the area fraction of the grid cell where the … … 228 247 !********************************************************** 229 248 249 species_loop: & 230 250 do ks=1,nspec ! loop over species 251 231 252 wetdeposit(ks)=0. 232 253 wetscav=0. 233 254 234 255 ! NIK 09.12.13: allowed to turn off either below cloud or in-cloud by including TWO separate blocks of if tests 235 236 237 !****i******************* 238 ! BELOW CLOUD SCAVENGING 239 !*********************** 240 if (weta(ks).gt.0. .and. clouds_v.ge.4) then !if positive below-cloud coefficient given in SPECIES file, and if in below cloud regime 241 ! for aerosols and not highliy soluble substances weta=5E-6 242 wetscav=weta(ks)*prec**wetb(ks) ! scavenging coefficient 243 endif !end below-cloud 244 245 246 !******************** 247 ! IN CLOUD SCAVENGING 248 !******************** 249 if (weta_in(ks).gt.0. .and. clouds_v.lt.4) then !if positive in-cloud coefficient given in SPECIES file, and if in in-cloud regime 250 251 if (ngrid.gt.0) then 252 act_temp=ttn(ix,jy,hz,n,ngrid) 253 else 254 act_temp=tt(ix,jy,hz,n) 255 endif 256 if (weta(ks).gt.0.) then 257 ! positive below-cloud coefficient from SPECIES file 258 259 if (indcloud .eq. 1) then ! below-cloud scavenging 260 261 ! for aerosols and not highly soluble substances weta=5E-6 262 wetscav=weta(ks)*prec**wetb(ks) ! scavenging coefficient 263 264 elseif (indcloud .eq. 2) then ! in-cloud scavenging 265 266 if (ngrid.gt.0) then 267 act_temp=ttn(ix,jy,kz,n,ngrid) 268 else 269 act_temp=tt(ix,jy,kz,n) 270 endif 256 271 257 272 ! NIK 31.01.2013: SPECIES defined parameters for the in-cloud scavening … … 259 274 ! wetb_in=0.36 (default) 260 275 ! wetc_in=0.9 (default) 261 ! wetd_in: Scaling factor for the total in-cloud scavenging (default 1.0 -no scaling)276 ! wetd_in: Scaling factor for the total in-cloud scavenging (default 1.0 - no scaling) 262 277 cl=weta_in(ks)*prec**wetb_in(ks) 263 278 if (dquer(ks).gt.0) then ! is particle 264 279 S_i=wetc_in(ks)/cl 265 280 else ! is gas 266 cle=(1 -cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl267 S_i=1 /cle281 cle=(1.-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl 282 S_i=1./cle 268 283 endif 269 wetscav=S_i*prec/3.6E6/clouds_h/wetd_in(ks) 270 ! write(*,*) 'in. wetscav:' 271 ! + ,wetscav,cle,cl,act_temp,prec,clouds_h 272 273 endif ! end in-cloud 274 275 284 wetscav=S_i*prec/3.6E6/wetd_in(ks)/(ictop-icbot) ! 3.6e6 converts mm/h to m/s 285 ! PS - prevent extrem values of wetscav: 286 wetscavold = 2.e-5*prec**0.8 287 wetscav = min( 0.1*wetscav, wetscavold ) ! here we sacle the above wetscav. This is now twice, also in wetd_in. 288 289 else ! PS: no cloud diagnosed, old scheme, 290 291 ! PS using with fixed a,b for simplicity, one may wish to change!! 292 wetscav = 2.e-5*prec**0.8 ! was before: 1.e-4*prec**0.62 293 294 endif ! end regime cases 276 295 277 !************************************************** 278 ! CALCULATE DEPOSITION 279 !************************************************** 280 ! if (wetscav.le.0) write (*,*) 'neg, or 0 wetscav!' 281 ! + ,wetscav,cle,cl,act_temp,prec,clouds_h,clouds_v 282 283 if (wetscav.gt.0.) then 284 wetdeposit(ks)=xmass1(jpart,ks)* & 285 (1.-exp(-wetscav*abs(ltsample)))*grfraction ! wet deposition 286 else ! if no scavenging 287 wetdeposit(ks)=0. 288 endif 289 290 296 ! calculate deposition 297 wetdeposit(ks)=xmass1(jpart,ks)* & 298 (1.-exp(-wetscav*abs(ltsample)))*grfraction 291 299 restmass = xmass1(jpart,ks)-wetdeposit(ks) 292 300 if (ioutputforeachrelease.eq.1) then … … 297 305 if (restmass .gt. smallnum) then 298 306 xmass1(jpart,ks)=restmass 299 ! depostatistic 300 ! wetdepo_sum(ks,kp)=wetdepo_sum(ks,kp)+wetdeposit(ks) 301 ! depostatistic 307 ! depostatistic: 308 !! wetdepo_sum(ks,kp)=wetdepo_sum(ks,kp)+wetdeposit(ks) 302 309 else 303 310 xmass1(jpart,ks)=0. … … 305 312 ! Correct deposited mass to the last time step when radioactive decay of 306 313 ! gridded deposited mass was calculated 307 if (decay(ks).gt.0.) then314 if (decay(ks).gt.0.) & 308 315 wetdeposit(ks)=wetdeposit(ks)*exp(abs(ldeltat)*decay(ks)) 309 endif 310 311 312 313 end do !all species 314 315 ! Sabine Eckhardt, June 2008 create deposition runs only for forward runs 316 ! Add the wet deposition to accumulated amount on output grid and nested output grid 317 !***************************************************************************** 316 317 else ! weta(k) <= 0 318 319 wetdeposit(ks)=0. 320 321 endif 322 323 enddo species_loop 324 325 ! Sabine Eckhardt, June 2008: write deposition only in forward runs 326 ! add the wet deposition from this step to accumulated amount 327 ! on output grid and nested output grid 318 328 319 329 if (ldirect.eq.1) then 320 call wetdepokernel(nclass(jpart),wetdeposit, real(xtra1(jpart)),&321 real(ytra1(jpart)),nage,kp)330 call wetdepokernel(nclass(jpart),wetdeposit, & 331 real(xtra1(jpart)),real(ytra1(jpart)),itage,nage,kp) 322 332 if (nested_output.eq.1) call wetdepokernel_nest(nclass(jpart), & 323 wetdeposit,real(xtra1(jpart)),real(ytra1(jpart)), nage,kp)324 endif 325 326 20 continue 327 end do ! all particles333 wetdeposit,real(xtra1(jpart)),real(ytra1(jpart)),itage,nage,kp) 334 endif 335 336 20 continue ! jump here for particles not to be treated 337 enddo particle_loop 328 338 329 339 end subroutine wetdepo -
branches/petra/src/wetdepokernel.f90
r4 r37 1 1 !********************************************************************** 2 ! Copyright 1998 ,1999,2000,2001,2002,2005,2007,2008,2009,2010*2 ! Copyright 1998-2015 * 3 3 ! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, * 4 4 ! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann * … … 20 20 !********************************************************************** 21 21 22 subroutine wetdepokernel(nunc,deposit,x,y, nage,kp)23 ! i i i ii22 subroutine wetdepokernel(nunc,deposit,x,y,itage,nage,kp) 23 ! i i i i i i 24 24 !***************************************************************************** 25 25 ! * 26 26 ! Attribution of the deposition from an individual particle to the * 27 ! deposition fields using a uniform kernel with bandwidths dxout and dyout.* 27 ! deposition fields using a uniform kernel with bandwidths dxout * 28 ! and dyout. * 28 29 ! * 29 30 ! Author: A. Stohl * 30 31 ! * 31 32 ! 26 December 1996 * 33 ! 34 ! PS, 2/2015: do not use kernel for itage < 3 h 35 ! same as for concentration in conccalc.f90 32 36 ! * 33 37 !***************************************************************************** … … 48 52 49 53 real :: x,y,deposit(maxspec),ddx,ddy,xl,yl,wx,wy,w 50 integer :: ix,jy,ixp,jyp,nunc,nage,ks,kp 54 integer :: ix,jy,ixp,jyp,nunc,nage,ks,kp,itage 51 55 52 56 xl=(x*dx+xoutshift)/dxout … … 74 78 75 79 76 ! Determine mass fractions for four grid points 77 !********************************************** 80 if (itage .lt. itagekernmin) then 81 ! no kernel, direct attribution to grid cell 82 83 do ks=1,nspec 84 if (ix.ge.0 .and. jy.ge.0 .and. & 85 ix.le.numxgrid-1 .and. jy.le.numygrid-1) & 86 wetgridunc(ix,jy,ks,kp,nunc,nage)= & 87 wetgridunc(ix,jy,ks,kp,nunc,nage)+deposit(ks) 88 enddo 89 90 else 91 ! Determine mass fractions for four grid points & distribute 78 92 79 do ks=1,nspec93 do ks=1,nspec 80 94 81 if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. & 82 (jy.le.numygrid-1)) then 83 w=wx*wy 84 wetgridunc(ix,jy,ks,kp,nunc,nage)= & 85 wetgridunc(ix,jy,ks,kp,nunc,nage)+deposit(ks)*w 95 if (ix.ge.0 .and. jy.ge.0 .and. & 96 ix.le.numxgrid-1 .and. jy.le.numygrid-1) then 97 w=wx*wy 98 wetgridunc(ix,jy,ks,kp,nunc,nage)= & 99 wetgridunc(ix,jy,ks,kp,nunc,nage)+deposit(ks)*w 100 endif 101 102 if ((ixp.ge.0).and.(jyp.ge.0).and.& 103 (ixp.le.numxgrid-1).and. (jyp.le.numygrid-1)) then 104 w=(1.-wx)*(1.-wy) 105 wetgridunc(ixp,jyp,ks,kp,nunc,nage)= & 106 wetgridunc(ixp,jyp,ks,kp,nunc,nage)+deposit(ks)*w 107 endif 108 109 if ((ixp.ge.0).and.(jy.ge.0).and. & 110 (ixp.le.numxgrid-1).and.(jy.le.numygrid-1)) then 111 w=(1.-wx)*wy 112 wetgridunc(ixp,jy,ks,kp,nunc,nage)= & 113 wetgridunc(ixp,jy,ks,kp,nunc,nage)+deposit(ks)*w 114 endif 115 116 if ((ix.ge.0).and.(jyp.ge.0).and. & 117 (ix.le.numxgrid-1).and.(jyp.le.numygrid-1)) then 118 w=wx*(1.-wy) 119 wetgridunc(ix,jyp,ks,kp,nunc,nage)= & 120 wetgridunc(ix,jyp,ks,kp,nunc,nage)+deposit(ks)*w 121 endif 122 123 end do 124 86 125 endif 87 126 88 if ((ixp.ge.0).and.(jyp.ge.0).and.(ixp.le.numxgrid-1).and. &89 (jyp.le.numygrid-1)) then90 w=(1.-wx)*(1.-wy)91 wetgridunc(ixp,jyp,ks,kp,nunc,nage)= &92 wetgridunc(ixp,jyp,ks,kp,nunc,nage)+deposit(ks)*w93 endif94 95 if ((ixp.ge.0).and.(jy.ge.0).and.(ixp.le.numxgrid-1).and. &96 (jy.le.numygrid-1)) then97 w=(1.-wx)*wy98 wetgridunc(ixp,jy,ks,kp,nunc,nage)= &99 wetgridunc(ixp,jy,ks,kp,nunc,nage)+deposit(ks)*w100 endif101 102 if ((ix.ge.0).and.(jyp.ge.0).and.(ix.le.numxgrid-1).and. &103 (jyp.le.numygrid-1)) then104 w=wx*(1.-wy)105 wetgridunc(ix,jyp,ks,kp,nunc,nage)= &106 wetgridunc(ix,jyp,ks,kp,nunc,nage)+deposit(ks)*w107 endif108 end do109 110 127 end subroutine wetdepokernel -
branches/petra/src/wetdepokernel_nest.f90
r4 r37 1 1 !********************************************************************** 2 ! Copyright 1998 ,1999,2000,2001,2002,2005,2007,2008,2009,2010*2 ! Copyright 1998-2015 * 3 3 ! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, * 4 4 ! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann * … … 20 20 !********************************************************************** 21 21 22 subroutine wetdepokernel_nest & 23 (nunc,deposit,x,y,nage,kp) 24 ! i i i i i i 22 subroutine wetdepokernel_nest (nunc,deposit,x,y,itage,nage,kp) 23 ! i i i i i i i 25 24 !***************************************************************************** 26 25 ! * … … 35 34 ! 2 September 2004: Adaptation from wetdepokernel. * 36 35 ! * 36 ! 37 ! PS, 2/2015: do not use kernel for itage < 3 h 38 ! same as for concentration in conccalc.f90 37 39 ! * 38 40 !***************************************************************************** … … 53 55 54 56 real :: x,y,deposit(maxspec),ddx,ddy,xl,yl,wx,wy,w 55 integer :: ix,jy,ixp,jyp,ks,kp,nunc,nage 57 integer :: ix,jy,ixp,jyp,ks,kp,nunc,nage,itage 56 58 57 59 … … 80 82 endif 81 83 84 if (itage .lt. itagekernmin) then 85 ! no kernel, direct attribution to grid cell 86 87 do ks=1,nspec 88 if (ix.ge.0 .and. jy.ge.0 .and. & 89 ix.le.numxgridn-1 .and. jy.le.numygridn-1) & 90 wetgriduncn(ix,jy,ks,kp,nunc,nage)= & 91 wetgriduncn(ix,jy,ks,kp,nunc,nage)+deposit(ks) 92 enddo 93 94 else 95 ! Determine mass fractions for four grid points & distribute 82 96 83 ! Determine mass fractions for four grid points 84 !********************************************** 97 do ks=1,nspec 85 98 86 do ks=1,nspec 99 if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgridn-1).and. & 100 (jy.le.numygridn-1)) then 101 w=wx*wy 102 wetgriduncn(ix,jy,ks,kp,nunc,nage)= & 103 wetgriduncn(ix,jy,ks,kp,nunc,nage)+deposit(ks)*w 104 endif 87 105 88 if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgridn-1).and. & 89 (jy.le.numygridn-1)) then 90 w=wx*wy 91 wetgriduncn(ix,jy,ks,kp,nunc,nage)= & 92 wetgriduncn(ix,jy,ks,kp,nunc,nage)+deposit(ks)*w 106 if ((ixp.ge.0).and.(jyp.ge.0).and.(ixp.le.numxgridn-1).and. & 107 (jyp.le.numygridn-1)) then 108 w=(1.-wx)*(1.-wy) 109 wetgriduncn(ixp,jyp,ks,kp,nunc,nage)= & 110 wetgriduncn(ixp,jyp,ks,kp,nunc,nage)+deposit(ks)*w 111 endif 112 113 if ((ixp.ge.0).and.(jy.ge.0).and.(ixp.le.numxgridn-1).and. & 114 (jy.le.numygridn-1)) then 115 w=(1.-wx)*wy 116 wetgriduncn(ixp,jy,ks,kp,nunc,nage)= & 117 wetgriduncn(ixp,jy,ks,kp,nunc,nage)+deposit(ks)*w 118 endif 119 120 if ((ix.ge.0).and.(jyp.ge.0).and.(ix.le.numxgridn-1).and. & 121 (jyp.le.numygridn-1)) then 122 w=wx*(1.-wy) 123 wetgriduncn(ix,jyp,ks,kp,nunc,nage)= & 124 wetgriduncn(ix,jyp,ks,kp,nunc,nage)+deposit(ks)*w 125 endif 126 127 end do 128 93 129 endif 94 95 if ((ixp.ge.0).and.(jyp.ge.0).and.(ixp.le.numxgridn-1).and. &96 (jyp.le.numygridn-1)) then97 w=(1.-wx)*(1.-wy)98 wetgriduncn(ixp,jyp,ks,kp,nunc,nage)= &99 wetgriduncn(ixp,jyp,ks,kp,nunc,nage)+deposit(ks)*w100 endif101 102 if ((ixp.ge.0).and.(jy.ge.0).and.(ixp.le.numxgridn-1).and. &103 (jy.le.numygridn-1)) then104 w=(1.-wx)*wy105 wetgriduncn(ixp,jy,ks,kp,nunc,nage)= &106 wetgriduncn(ixp,jy,ks,kp,nunc,nage)+deposit(ks)*w107 endif108 109 if ((ix.ge.0).and.(jyp.ge.0).and.(ix.le.numxgridn-1).and. &110 (jyp.le.numygridn-1)) then111 w=wx*(1.-wy)112 wetgriduncn(ix,jyp,ks,kp,nunc,nage)= &113 wetgriduncn(ix,jyp,ks,kp,nunc,nage)+deposit(ks)*w114 endif115 116 end do117 130 end subroutine wetdepokernel_nest
Note: See TracChangeset
for help on using the changeset viewer.