Changeset f251e57 in flexpart.git
- Timestamp:
- Jun 25, 2018, 7:44:14 PM (6 years ago)
- Branches:
- univie
- Children:
- c0884a8
- Parents:
- 8b3d324
- Location:
- src
- Files:
-
- 3 edited
-
changelog.txt (modified) (2 diffs)
-
verttransform_ecmwf.f90 (modified) (4 diffs)
-
verttransform_gfs.f90 (modified) (23 diffs)
Legend:
- Unmodified
- Added
- Removed
-
src/changelog.txt
r1a8fbee rf251e57 50 50 7) Bring some more code under the IF (INIT) block 51 51 8) make some annotations 52 [minor modification PS 2018-06-25] 52 53 53 54 * com_mod.f90 PS 2018-06-13 … … 63 64 3) Change name of loops to represent the index 64 65 4) make some annotations 66 67 * verttransform_ecmwf.f90 PS 2018-06-25 68 =================== 69 70 1) Some code layout improvement in the first part 71 2) Fix ticket:140 Introduce a new way of establishing the reference position 72 for the vertical grid. Also correct a minor bug in the calculation 73 of height (array assignment where it was not intended) 74 3) Add back the SAVE attribute to INIT, just to be sure -
src/verttransform_ecmwf.f90
r8b3d324 rf251e57 157 157 158 158 if (ioldref .eq. 0) then ! old reference grid point 159 160 161 162 163 164 165 166 167 159 do jy=0,nymin1 160 do ix=0,nxmin1 161 if (ps(ix,jy,1,n).gt.1000.e2) then 162 ixref=ix 163 jyref=jy 164 goto 3 165 endif 166 end do 167 end do 168 168 3 continue 169 169 … … 175 175 176 176 psmean = sum( ps(:,:,1,n) ) / (nx*ny) 177 print*,'height: fg psmean',psmean177 ! print*,'height: fg psmean',psmean 178 178 psstd = sqrt(sum( (ps(:,:,1,n) - psmean)**2 ) / (nx*ny)) 179 179 … … 201 201 202 202 pintref = akz(kz)+bkz(kz)*ps(ixref,jyref,1,n) 203 tv = tth(ixref,jyref,kz,n)*(1.+0.608*qvh(ixref,jyref,kz,n))203 tv(ixref,jyref) = tth(ixref,jyref,kz,n)*(1.+0.608*qvh(ixref,jyref,kz,n)) 204 204 205 205 if (abs(tv(ixref,jyref)-tvoldref) .gt. 0.2) then … … 214 214 tvoldref = tv(ixref,jyref) 215 215 poldref = pintref 216 print*,'height=',kz,height(kz),tvoldref,poldref216 ! print*,'height=',kz,height(kz),tvoldref,poldref 217 217 218 218 end do -
src/verttransform_gfs.f90
r6ecb30a rf251e57 21 21 22 22 subroutine verttransform_gfs(n,uuh,vvh,wwh,pvh) 23 ! i i i i i 24 !***************************************************************************** 25 ! * 26 ! This subroutine transforms temperature, dew point temperature and * 27 ! wind components from eta to meter coordinates. * 28 ! The vertical wind component is transformed from Pa/s to m/s using * 29 ! the conversion factor pinmconv. * 30 ! In addition, this routine calculates vertical density gradients * 31 ! needed for the parameterization of the turbulent velocities. * 32 ! * 33 ! Author: A. Stohl, G. Wotawa * 34 ! * 35 ! 12 August 1996 * 36 ! Update: 16 January 1998 * 37 ! * 38 ! Major update: 17 February 1999 * 39 ! by G. Wotawa * 40 ! CHANGE 17/11/2005 Caroline Forster, NCEP GFS version * 41 ! * 42 ! - Vertical levels for u, v and w are put together * 43 ! - Slope correction for vertical velocity: Modification of calculation * 44 ! procedure * 45 ! * 46 !***************************************************************************** 47 ! Changes, Bernd C. Krueger, Feb. 2001: 48 ! Variables tth and qvh (on eta coordinates) from common block 49 ! 50 ! Unified ECMWF and GFS builds 51 ! Marian Harustak, 12.5.2017 52 ! - Renamed routine from verttransform to verttransform_gfs 53 ! 54 !***************************************************************************** 55 ! * 56 ! Variables: * 57 ! nx,ny,nz field dimensions in x,y and z direction * 58 ! uu(0:nxmax,0:nymax,nzmax,2) wind components in x-direction [m/s] * 59 ! vv(0:nxmax,0:nymax,nzmax,2) wind components in y-direction [m/s] * 60 ! ww(0:nxmax,0:nymax,nzmax,2) wind components in z-direction [deltaeta/s]* 61 ! tt(0:nxmax,0:nymax,nzmax,2) temperature [K] * 62 ! pv(0:nxmax,0:nymax,nzmax,2) potential voriticity (pvu) * 63 ! ps(0:nxmax,0:nymax,2) surface pressure [Pa] * 64 ! clouds(0:nxmax,0:nymax,0:nzmax,2) cloud field for wet deposition * 65 ! * 66 !***************************************************************************** 23 ! i i i i i 24 !***************************************************************************** 25 ! * 26 ! This subroutine transforms temperature, dew point temperature and * 27 ! wind components from eta to meter coordinates. * 28 ! The vertical wind component is transformed from Pa/s to m/s using * 29 ! the conversion factor pinmconv. * 30 ! In addition, this routine calculates vertical density gradients * 31 ! needed for the parameterization of the turbulent velocities. * 32 ! * 33 ! Author: A. Stohl, G. Wotawa * 34 ! * 35 ! 12 August 1996 * 36 ! Update: 16 January 1998 * 37 ! * 38 ! Major update: 17 February 1999 * 39 ! by G. Wotawa * 40 ! CHANGE 17/11/2005 Caroline Forster, NCEP GFS version * 41 ! * 42 ! - Vertical levels for u, v and w are put together * 43 ! - Slope correction for vertical velocity: Modification of calculation * 44 ! procedure * 45 ! * 46 !***************************************************************************** 47 ! Changes, Bernd C. Krueger, Feb. 2001: 48 ! Variables tth and qvh (on eta coordinates) from common block 49 ! 50 ! Unified ECMWF and GFS builds 51 ! Marian Harustak, 12.5.2017 52 ! - Renamed from verttransform to verttransform_ecmwf 53 ! 54 ! undocumented modifications by NILU for v10 * 55 ! * 56 ! Petra Seibert, 2018-06-13: * 57 ! - fix bug of ticket:140 (find reference position for vertical grid) * 58 ! - put back SAVE attribute for INIT, just to be safe * 59 ! - minor changes, most of them just cosmetics * 60 ! for details see changelog.txt * 61 ! * 62 !***************************************************************************** 63 ! * 64 ! Variables: * 65 ! nx,ny,nz field dimensions in x,y and z direction * 66 ! uu(0:nxmax,0:nymax,nzmax,2) wind components in x-direction [m/s] * 67 ! vv(0:nxmax,0:nymax,nzmax,2) wind components in y-direction [m/s] * 68 ! ww(0:nxmax,0:nymax,nzmax,2) wind components in z-direction [deltaeta/s]* 69 ! tt(0:nxmax,0:nymax,nzmax,2) temperature [K] * 70 ! pv(0:nxmax,0:nymax,nzmax,2) potential voriticity (pvu) * 71 ! ps(0:nxmax,0:nymax,2) surface pressure [Pa] * 72 ! clouds(0:nxmax,0:nymax,0:nzmax,2) cloud field for wet deposition * 73 ! * 74 !***************************************************************************** 67 75 68 76 use par_mod … … 72 80 implicit none 73 81 74 integer :: ix,jy,kz,iz,n,kmin,kl,klp,ix1,jy1,ixp,jyp ,ixm,jym82 integer :: ix,jy,kz,iz,n,kmin,kl,klp,ix1,jy1,ixp,jyp 75 83 integer :: rain_cloud_above,kz_inv 76 84 real :: f_qvsat,pressure … … 78 86 real :: rhoh(nuvzmax),pinmconv(nzmax) 79 87 real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi 88 89 !> for reference profile (PS) 90 real :: tvoldref, poldref, pintref, psmean, psstd 91 integer :: ixyref(2), ixref,jyref 92 integer, parameter :: ioldref = 1 ! PS 2018-06-13: set to 0 if you 93 !! want old method of searching reference location for the vertical grid 94 !! 1 for new method (options for other methods 2,... in the future) 95 80 96 real :: xlon,ylat,xlonr,dzdx,dzdy 81 97 real :: dzdx1,dzdx2,dzdy1,dzdy2,cosf … … 88 104 real,parameter :: const=r_air/ga 89 105 90 ! NCEP version 106 logical, save :: init = .true. ! PS 2018-06-13: add back save attribute 107 108 !> GFS version 91 109 integer :: llev, i 92 110 93 logical :: init = .true. 94 95 96 !************************************************************************* 97 ! If verttransform is called the first time, initialize heights of the * 98 ! z levels in meter. The heights are the heights of model levels, where * 99 ! u,v,T and qv are given. * 100 !************************************************************************* 111 112 113 !************************************************************************* 114 ! If verttransform is called the first time, initialize heights of the * 115 ! z levels in meter. The heights are the heights of model levels, where * 116 ! u,v,T and qv are given. * 117 !************************************************************************* 101 118 102 119 if (init) then 103 120 104 ! Search for a point with high surface pressure (i.e. not above significant topography) 105 ! Then, use this point to construct a reference z profile, to be used at all times 106 !***************************************************************************** 107 108 do jy=0,nymin1 109 do ix=0,nxmin1 110 if (ps(ix,jy,1,n).gt.100000.) then 111 ixm=ix 112 jym=jy 121 ! Search for a point with high surface pressure (i.e. not above significant topography) 122 ! Then, use this point to construct a reference z profile, to be used at all times 123 !***************************************************************************** 124 125 if (ioldref .eq. 0) then ! old reference grid point 126 do jy=0,nymin1 127 do ix=0,nxmin1 128 if (ps(ix,jy,1,n).gt.1000.e2) then 129 ixref=ix 130 jyref=jy 113 131 goto 3 114 132 endif … … 116 134 end do 117 135 3 continue 118 119 120 tvold=tt2(ixm,jym,1,n)*(1.+0.378*ew(td2(ixm,jym,1,n))/ & 121 ps(ixm,jym,1,n)) 122 pold=ps(ixm,jym,1,n) 136 ! print*,'oldheights at' ,ixref,jyref,ps(ixref,jyref,1,n) 137 else ! new reference grid point 138 ! PS: the old version fails if the pressure is <=1000 hPa in the whole 139 ! domain. Let us find a good replacement, not just a quick fix. 140 ! Search point near to mean pressure after excluding mountains 141 142 psmean = sum( ps(:,:,1,n) ) / (nx*ny) 143 ! print*,'height: fg psmean',psmean 144 psstd = sqrt(sum( (ps(:,:,1,n) - psmean)**2 ) / (nx*ny)) 145 146 !> new psmean using only points within $\plusminus\sigma$ 147 ! psmean = sum( ps(:,:,1,n), abs(ps(:,:,1,n)-psmean) < psstd ) / & 148 ! count(abs(ps(:,:,1,n)-psmean) < psstd) 149 150 !> new psmean using only points with $p\gt \overbar{p}+\sigma_p$ 151 !> (reject mountains, accept valleys) 152 psmean = sum( ps(:,:,1,n), ps(:,:,1,n) > psmean - psstd ) / & 153 count(ps(:,:,1,n) > psmean - psstd) 154 ! print*,'height: std, new psmean',psstd,psmean 155 ixyref = minloc( abs( ps(:,:,1,n) - psmean ) ) 156 ixref = ixyref(1) 157 jyref = ixyref(2) 158 ! print*,'newheights at' ,ixref,jyref,ps(ixref,jyref,1,n) 159 endif 160 161 tvoldref=tt2(ixref,jyref,1,n)* & 162 ( 1. + 0.378*ew(td2(ixref,jyref,1,n) ) / ps(ixref,jyref,1,n)) 163 poldref=ps(ixref,jyref,1,n) 123 164 height(1)=0. 165 kz=1 166 ! print*,'height=',kz,height(kz),tvoldref,poldref 124 167 125 168 do kz=2,nuvz 126 pint=akz(kz)+bkz(kz)*ps(ixm,jym,1,n) 127 tv=tth(ixm,jym,kz,n)*(1.+0.608*qvh(ixm,jym,kz,n)) 128 129 if (abs(tv-tvold).gt.0.2) then 130 height(kz)=height(kz-1)+const*log(pold/pint)* & 131 (tv-tvold)/log(tv/tvold) 169 170 pintref = akz(kz) 171 ! Note that for GFS data, the akz variable contains the input level 172 ! pressure values. bkz is zero (I removed all terms with bkz that 173 ! were erroneously copied from verttransform_ecmwf). [PS, June 2018] 174 175 tv = tth(ixref,jyref,kz,n)*(1.+0.608*qvh(ixref,jyref,kz,n)) 176 177 if (abs(tv-tvoldref) .gt. 0.2) then 178 height(kz)=height(kz-1) + & 179 const * log(poldref/pintref) * (tv-tvoldref) / log(tv/tvoldref) 132 180 else 133 height(kz) =height(kz-1)+const*log(pold/pint)*tv181 height(kz) = height(kz-1) + const*log(poldref/pintref)*tv 134 182 endif 135 183 136 tvold=tv 137 pold=pint 138 end do 139 140 141 ! Determine highest levels that can be within PBL 142 !************************************************ 184 tvoldref=tv 185 poldref=pintref 186 ! print*,'height=',kz,height(kz),tvoldref,poldref 187 end do 188 189 190 ! Determine highest levels that can be within PBL 191 !************************************************ 143 192 144 193 do kz=1,nz … … 150 199 9 continue 151 200 152 153 201 ! Do not repeat initialization of the Cartesian z grid 202 !***************************************************** 154 203 155 204 init=.false. 156 205 157 endif 158 159 160 161 206 endif ! init block (vertical grid construction) 207 208 209 ! Loop over the whole grid 210 !************************* 162 211 163 212 do jy=0,nymin1 164 213 do ix=0,nxmin1 165 214 166 215 ! NCEP version: find first level above ground 167 216 llev = 0 168 217 do i=1,nuvz … … 171 220 llev = llev+1 172 221 if (llev.gt.nuvz-2) llev = nuvz-2 173 174 175 176 177 178 179 222 ! if (llev.eq.nuvz-2) write(*,*) 'verttransform 223 ! +WARNING: LLEV eq NUZV-2' 224 ! NCEP version 225 226 227 ! compute height of pressure levels above ground 228 !*********************************************** 180 229 181 230 tvold=tth(ix,jy,llev,n)*(1.+0.608*qvh(ix,jy,llev,n)) … … 202 251 end do 203 252 204 253 ! pinmconv=(h2-h1)/(p2-p1) 205 254 206 255 pinmconv(llev)=(uvwzlev(ix,jy,llev+1)-uvwzlev(ix,jy,llev))/ & … … 217 266 218 267 219 220 268 ! Levels, where u,v,t and q are given 269 !************************************ 221 270 222 271 uu(ix,jy,1,n)=uuh(ix,jy,llev) … … 267 316 268 317 269 270 318 ! Levels, where w is given 319 !************************* 271 320 272 321 ww(ix,jy,1,n)=wwh(ix,jy,llev)*pinmconv(llev) … … 287 336 288 337 289 290 338 ! Compute density gradients at intermediate levels 339 !************************************************* 291 340 292 341 drhodz(ix,jy,1,n)=(rho(ix,jy,2,n)-rho(ix,jy,1,n))/ & … … 302 351 303 352 304 305 306 307 353 !**************************************************************** 354 ! Compute slope of eta levels in windward direction and resulting 355 ! vertical wind correction 356 !**************************************************************** 308 357 309 358 do jy=1,ny-2 … … 311 360 do ix=1,nx-2 312 361 313 362 ! NCEP version: find first level above ground 314 363 llev = 0 315 364 do i=1,nuvz … … 318 367 llev = llev+1 319 368 if (llev.gt.nuvz-2) llev = nuvz-2 320 321 322 369 ! if (llev.eq.nuvz-2) write(*,*) 'verttransform 370 ! +WARNING: LLEV eq NUZV-2' 371 ! NCEP version 323 372 324 373 kmin=llev+1 … … 361 410 362 411 363 364 365 412 ! If north pole is in the domain, calculate wind velocities in polar 413 ! stereographic coordinates 414 !******************************************************************* 366 415 367 416 if (nglobal) then … … 381 430 do iz=1,nz 382 431 383 432 ! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT 384 433 xlon=xlon0+real(nx/2-1)*dx 385 434 xlonr=xlon*pi/180. … … 396 445 if(ddpol.gt.2.0*pi) ddpol=ddpol-2.0*pi 397 446 398 447 ! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID 399 448 xlon=180.0 400 449 xlonr=xlon*pi/180. … … 411 460 412 461 413 414 462 ! Fix: Set W at pole to the zonally averaged W of the next equator- 463 ! ward parallel of latitude 415 464 416 465 do iz=1,nz … … 430 479 431 480 432 433 434 481 ! If south pole is in the domain, calculate wind velocities in polar 482 ! stereographic coordinates 483 !******************************************************************* 435 484 436 485 if (sglobal) then … … 448 497 do iz=1,nz 449 498 450 499 ! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT 451 500 xlon=xlon0+real(nx/2-1)*dx 452 501 xlonr=xlon*pi/180. … … 462 511 if(ddpol.gt.2.0*pi) ddpol=ddpol-2.0*pi 463 512 464 513 ! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID 465 514 xlon=180.0 466 515 xlonr=xlon*pi/180. … … 478 527 479 528 480 481 529 ! Fix: Set W at pole to the zonally averaged W of the next equator- 530 ! ward parallel of latitude 482 531 483 532 do iz=1,nz … … 496 545 497 546 498 499 500 547 ! write (*,*) 'initializing clouds, n:',n,nymin1,nxmin1,nz 548 ! create a cloud and rainout/washout field, clouds occur where rh>80% 549 ! total cloudheight is stored at level 0 501 550 do jy=0,nymin1 502 551 do ix=0,nxmin1
Note: See TracChangeset
for help on using the changeset viewer.
