[16] | 1 | !*********************************************************************** |
---|
| 2 | !* Copyright 2012,2013 * |
---|
| 3 | !* Jerome Brioude, Delia Arnold, Andreas Stohl, Wayne Angevine, * |
---|
| 4 | !* John Burkhart, Massimo Cassiani, Adam Dingwell, Richard C Easter, Sabine Eckhardt,* |
---|
| 5 | !* Stephanie Evan, Jerome D Fast, Don Morton, Ignacio Pisso, * |
---|
| 6 | !* Petra Seibert, Gerard Wotawa, Caroline Forster, Harald Sodemann, * |
---|
| 7 | !* * |
---|
| 8 | !* This file is part of FLEXPART WRF * |
---|
| 9 | !* * |
---|
| 10 | !* FLEXPART is free software: you can redistribute it and/or modify * |
---|
| 11 | !* it under the terms of the GNU General Public License as published by* |
---|
| 12 | !* the Free Software Foundation, either version 3 of the License, or * |
---|
| 13 | !* (at your option) any later version. * |
---|
| 14 | !* * |
---|
| 15 | !* FLEXPART is distributed in the hope that it will be useful, * |
---|
| 16 | !* but WITHOUT ANY WARRANTY; without even the implied warranty of * |
---|
| 17 | !* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * |
---|
| 18 | !* GNU General Public License for more details. * |
---|
| 19 | !* * |
---|
| 20 | !* You should have received a copy of the GNU General Public License * |
---|
| 21 | !* along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. * |
---|
| 22 | !*********************************************************************** |
---|
| 23 | |
---|
| 24 | subroutine readwind(indj,n,uuh,vvh,wwh,divh) |
---|
| 25 | !********************************************************************** |
---|
| 26 | ! * |
---|
| 27 | ! TRAJECTORY MODEL SUBROUTINE READWIND * |
---|
| 28 | ! * |
---|
| 29 | !********************************************************************** |
---|
| 30 | ! * |
---|
| 31 | ! AUTHOR: G. WOTAWA * |
---|
| 32 | ! DATE: 1997-08-05 * |
---|
| 33 | ! LAST UPDATE: 2000-10-17, Andreas Stohl * |
---|
| 34 | ! * |
---|
| 35 | ! Bernd C. Krueger, Feb. 2001: Variables tth and qvh * |
---|
| 36 | ! (on eta coordinates) in common block * |
---|
| 37 | ! * |
---|
| 38 | ! Oct-Dec, 2005: R. Easter. Major changes for WRF. * |
---|
| 39 | ! 06-nov-2005 rce - change uuh,vvh dimension back to original * |
---|
| 40 | ! 16-nov-2005 rce - zzh is shifted like pph,tth * |
---|
| 41 | ! * |
---|
| 42 | ! 11-June-2007, W.WANG -- read TKE, change ndims_exp=1 for P_TOP |
---|
| 43 | ! 19-Oct -2007, read tcc, RAINC, RAINNC, CLDFRA,W0AVG |
---|
| 44 | ! Note RAINC, RAINNC are accumulated prec |
---|
| 45 | ! Dec 2011, J Brioude: modifications, notably for the mean wind |
---|
| 46 | ! |
---|
| 47 | ! D. Arnold May 2012: quick fix - for CLDFRA |
---|
| 48 | ! CLDFRA was for sub-grid variability of precipitation |
---|
| 49 | ! within a cell from Hertel et al., 1995 (grid cells of approx. 150 km) |
---|
| 50 | ! WRF will be typically be run at higher resolutions, therefore |
---|
| 51 | ! we simply skip this condition and consider the maximum fraction |
---|
| 52 | ! no more reading of CLDFRA needed, initialized to 0. |
---|
| 53 | ! Same modifications for the nested domains! |
---|
| 54 | |
---|
| 55 | !********************************************************************** |
---|
| 56 | ! * |
---|
| 57 | ! Note: This is the FLEXPART_WRF version of subroutine readwind. * |
---|
| 58 | ! The met fields are read from WRF netcdf output files. * |
---|
| 59 | ! There are many differences from the FLEXPART version. * |
---|
| 60 | ! * |
---|
| 61 | ! DESCRIPTION: * |
---|
| 62 | ! * |
---|
| 63 | ! READING OF ECMWF METEOROLOGICAL FIELDS FROM INPUT DATA FILES. THE * |
---|
| 64 | ! INPUT DATA FILES ARE EXPECTED TO BE AVAILABLE IN GRIB CODE * |
---|
| 65 | ! * |
---|
| 66 | ! INPUT: * |
---|
| 67 | ! indj indicates number of the wind field to be read in * |
---|
| 68 | ! n temporal index for meteorological fields (1 to 3)* |
---|
| 69 | ! * |
---|
| 70 | ! IMPORTANT VARIABLES FROM COMMON BLOCK: * |
---|
| 71 | ! * |
---|
| 72 | ! wfname File name of data to be read in * |
---|
| 73 | ! nx,ny,nuvz,nwz expected field dimensions * |
---|
| 74 | ! nlev_ec number of "T-grid" vertical levels wwf model * |
---|
| 75 | ! (the unstaggered "bottom_top" dimension) * |
---|
| 76 | ! uu,vv,ww wind fields * |
---|
| 77 | ! tt,qv temperature and specific humidity * |
---|
| 78 | ! ps surface pressure * |
---|
| 79 | ! * |
---|
| 80 | !********************************************************************** |
---|
| 81 | ! |
---|
| 82 | |
---|
| 83 | ! include 'includepar' |
---|
| 84 | ! include 'includecom' |
---|
| 85 | use par_mod |
---|
| 86 | use com_mod |
---|
| 87 | |
---|
| 88 | ! subr arguments |
---|
| 89 | integer :: indj, n |
---|
| 90 | |
---|
| 91 | real(kind=4) :: uuh(0:nxmax-1,0:nymax-1,nuvzmax) |
---|
| 92 | real(kind=4) :: vvh(0:nxmax-1,0:nymax-1,nuvzmax) |
---|
| 93 | real(kind=4) :: wwh(0:nxmax-1,0:nymax-1,nwzmax) |
---|
| 94 | ! real(kind=4) :: urot(0:nxmax-1,0:nymax-1,nuvzmax) |
---|
| 95 | ! real(kind=4) :: vrot(0:nxmax-1,0:nymax-1,nuvzmax) |
---|
| 96 | real(kind=4) :: divh(0:nxmax-1,0:nymax-1,nuvzmax) |
---|
| 97 | real(kind=4) :: mu(0:nxmax-1,0:nymax-1,1),mu2 |
---|
| 98 | real(kind=4) :: mub(0:nxmax-1,0:nymax-1,1) |
---|
| 99 | ! real :: utrue1,vtrue1,utrue2,vtrue2,dumy |
---|
| 100 | ! real(kind=4) :: m_u(0:nxmax-1,0:nymax-1,1) |
---|
| 101 | ! real(kind=4) :: m_v(0:nxmax-1,0:nymax-1,1) |
---|
| 102 | ! real(kind=4) :: m_w(0:nxmax-1,0:nymax-1,1) |
---|
| 103 | |
---|
| 104 | ! real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax) |
---|
| 105 | ! real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax) |
---|
| 106 | ! real :: urot(0:nxmax-1,0:nymax-1,nuvzmax) |
---|
| 107 | ! real :: vrot(0:nxmax-1,0:nymax-1,nuvzmax) |
---|
| 108 | ! real :: wwh(0:nxmax-1,0:nymax-1,nwzmax) |
---|
| 109 | ! real :: divh(0:nxmax-1,0:nymax-1,nuvzmax) |
---|
| 110 | ! real :: mu(0:nxmax-1,0:nymax-1,1),mu2 |
---|
| 111 | ! real :: mub(0:nxmax-1,0:nymax-1,1) |
---|
| 112 | ! real :: m_u(0:nxmax-1,0:nymax-1,1) |
---|
| 113 | ! real :: m_v(0:nxmax-1,0:nymax-1,1) |
---|
| 114 | ! real :: m_x(0:nxmax-1,0:nymax-1,1) |
---|
| 115 | ! real :: m_y(0:nxmax-1,0:nymax-1,1) |
---|
| 116 | ! real :: m_z(0:nxmax-1,0:nymax-1,1) |
---|
| 117 | |
---|
| 118 | ! local variables |
---|
| 119 | integer,parameter :: ndims_max=4 |
---|
| 120 | |
---|
| 121 | integer :: i, idiagaa, ierr, ifn, itime |
---|
| 122 | integer :: iduma |
---|
| 123 | integer :: j, jhhmmss, jyyyymmdd |
---|
| 124 | integer :: k, kbgn |
---|
| 125 | integer :: lendim(ndims_max), lendim_exp(ndims_max), & |
---|
| 126 | lendim_max(ndims_max) |
---|
| 127 | integer :: levdiff2 |
---|
| 128 | integer :: ndims, ndims_exp |
---|
| 129 | integer :: n_west_east, n_south_north, n_bottom_top |
---|
| 130 | integer :: m_grid_id_dum, m_parent_grid_id_dum, & |
---|
| 131 | m_parent_grid_ratio_dum, & |
---|
| 132 | i_parent_start_dum, j_parent_start_dum, & |
---|
| 133 | map_proj_id_dum, & |
---|
| 134 | ext_scalar,pbl_physics,mp_physics_dum |
---|
| 135 | |
---|
| 136 | real :: dx_met, dy_met |
---|
| 137 | real :: duma, dumb, dumc, dumd, dume |
---|
| 138 | real :: dumdz |
---|
| 139 | ! real(kind=4) :: dumarray_aa(nwzmax+1) |
---|
| 140 | ! real(kind=4) :: dumarray_pp(0:nxmax-1,0:nymax-1,nwzmax+1) |
---|
| 141 | real :: dumarray_aa(nwzmax+1) |
---|
| 142 | real :: dumarray_pp(0:nxmax-1,0:nymax-1,nwzmax+1) |
---|
| 143 | real :: ewater_mb, esatwater_mb |
---|
| 144 | real :: ew ! this is an external function |
---|
| 145 | real :: map_stdlon_dum, map_truelat1_dum, map_truelat2_dum |
---|
| 146 | real :: pint |
---|
| 147 | real :: toler |
---|
| 148 | |
---|
| 149 | ! real(kind=4) :: ewss(0:nxmax-1,0:nymax-1),nsss(0:nxmax-1,0:nymax-1) |
---|
| 150 | real :: ewss(0:nxmax-1,0:nymax-1),nsss(0:nxmax-1,0:nymax-1) |
---|
| 151 | real :: plev1,pmean,tv,fu,hlev1,ff10m,fflev1 |
---|
| 152 | |
---|
| 153 | real(kind=dp) :: jul,juldate |
---|
| 154 | |
---|
| 155 | character(len=160) :: fnamenc, varname,fnamenc2 |
---|
| 156 | |
---|
| 157 | logical :: hflswitch |
---|
| 158 | |
---|
| 159 | ! |
---|
| 160 | ! get grid info from the wrf netcdf file |
---|
| 161 | ! and check it for consistency against values from gridcheck |
---|
| 162 | ! |
---|
| 163 | fnamenc = path(2)(1:length(2))//wfname(indj) |
---|
| 164 | idiagaa = 0 |
---|
| 165 | |
---|
| 166 | call read_ncwrfout_gridinfo( ierr, idiagaa, fnamenc, & |
---|
| 167 | n_west_east, n_south_north, n_bottom_top, & |
---|
| 168 | dx_met, dy_met, & |
---|
| 169 | m_grid_id_dum, m_parent_grid_id_dum, m_parent_grid_ratio_dum, & |
---|
| 170 | i_parent_start_dum, j_parent_start_dum, & |
---|
| 171 | map_proj_id_dum, map_stdlon_dum, & |
---|
| 172 | map_truelat1_dum, map_truelat2_dum, & |
---|
| 173 | ext_scalar,pbl_physics,mp_physics_dum ) |
---|
| 174 | if (ierr .ne. 0) then |
---|
| 175 | write(*,9100) 'error getting gridinfor for met file', fnamenc |
---|
| 176 | stop |
---|
| 177 | end if |
---|
| 178 | |
---|
| 179 | 9100 format( / '*** readwind -- ', a ) |
---|
| 180 | 9110 format( / '*** readwind -- ', a, 1x, i8 / & |
---|
| 181 | 'file = ', a ) |
---|
| 182 | 9120 format( / '*** readwind -- ', a, 2(1x,i8) / & |
---|
| 183 | 'file = ', a ) |
---|
| 184 | 9130 format( / '*** readwind -- ', a, 3(1x,i8) / & |
---|
| 185 | 'file = ', a ) |
---|
| 186 | 9115 format( / '*** readwind -- ', a / a, 1x, i8 / & |
---|
| 187 | 'file = ', a ) |
---|
| 188 | 9125 format( / '*** readwind -- ', a / a, 2(1x,i8) / & |
---|
| 189 | 'file = ', a ) |
---|
| 190 | 9135 format( / '*** readwind -- ', a / a, 3(1x,i8) / & |
---|
| 191 | 'file = ', a ) |
---|
| 192 | |
---|
| 193 | toler = 2.0e-7 |
---|
| 194 | |
---|
| 195 | if (nx .ne. n_west_east) then |
---|
| 196 | write(*,9100) 'nx not consistent', fnamenc |
---|
| 197 | stop |
---|
| 198 | end if |
---|
| 199 | if (ny .ne. n_south_north) then |
---|
| 200 | write(*,9100) 'ny not consistent', fnamenc |
---|
| 201 | stop |
---|
| 202 | end if |
---|
| 203 | if (nlev_ec .ne. n_bottom_top) then |
---|
| 204 | write(*,9100) 'nlev_ec not consistent', fnamenc |
---|
| 205 | stop |
---|
| 206 | end if |
---|
| 207 | if (nwz .ne. n_bottom_top+1) then |
---|
| 208 | write(*,9100) 'nwz not consistent', fnamenc |
---|
| 209 | stop |
---|
| 210 | end if |
---|
| 211 | ! if (nuvz .ne. n_bottom_top+1) then |
---|
| 212 | ! write(*,9100) 'nuvz not consistent', fnamenc |
---|
| 213 | ! stop |
---|
| 214 | ! end if |
---|
| 215 | |
---|
| 216 | if (m_grid_id(0) .ne. m_grid_id_dum) then |
---|
| 217 | write(*,9100) 'm_grid_id not consistent', fnamenc |
---|
| 218 | write(*,*) m_grid_id(0), m_grid_id_dum |
---|
| 219 | stop |
---|
| 220 | end if |
---|
| 221 | if (m_parent_grid_id(0) .ne. m_parent_grid_id_dum) then |
---|
| 222 | write(*,9100) 'm_parent_grid_id not consistent', fnamenc |
---|
| 223 | stop |
---|
| 224 | end if |
---|
| 225 | if (m_parent_grid_ratio(0) .ne. m_parent_grid_ratio_dum) then |
---|
| 226 | write(*,9100) 'm_parent_grid_ratio not consistent', fnamenc |
---|
| 227 | stop |
---|
| 228 | end if |
---|
| 229 | if (i_parent_start(0) .ne. i_parent_start_dum) then |
---|
| 230 | write(*,9100) 'i_parent_start not consistent', fnamenc |
---|
| 231 | stop |
---|
| 232 | end if |
---|
| 233 | if (j_parent_start(0) .ne. j_parent_start_dum) then |
---|
| 234 | write(*,9100) 'j_parent_start not consistent', fnamenc |
---|
| 235 | stop |
---|
| 236 | end if |
---|
| 237 | |
---|
| 238 | if (abs(dx - dx_met) .gt. toler*abs(dx)) then |
---|
| 239 | write(*,9100) 'dx not consistent', fnamenc |
---|
| 240 | stop |
---|
| 241 | end if |
---|
| 242 | if (abs(dy - dy_met) .gt. toler*abs(dy)) then |
---|
| 243 | write(*,9100) 'dy not consistent', fnamenc |
---|
| 244 | stop |
---|
| 245 | end if |
---|
| 246 | |
---|
| 247 | ! locate the date/time in the file |
---|
| 248 | itime = 0 |
---|
| 249 | 1100 itime = itime + 1 |
---|
| 250 | call read_ncwrfout_1datetime( ierr, fnamenc, & |
---|
| 251 | itime, jyyyymmdd, jhhmmss ) |
---|
| 252 | if (ierr .eq. -1) then |
---|
| 253 | write(*,9100) 'error reading time from met file', fnamenc |
---|
| 254 | stop |
---|
| 255 | else if (ierr .ne. 0) then |
---|
| 256 | write(*,9125) 'unable to locate date/time in met file', & |
---|
| 257 | 'indj, itime =', indj, itime, fnamenc |
---|
| 258 | stop |
---|
| 259 | else |
---|
| 260 | jul = juldate( jyyyymmdd, jhhmmss ) |
---|
| 261 | duma = (jul-bdate)*86400. |
---|
| 262 | iduma = nint(duma) |
---|
| 263 | if (iduma .ne. wftime(indj)) goto 1100 |
---|
| 264 | end if |
---|
| 265 | if (option_verbose.eq.1) then |
---|
| 266 | write(*,*) |
---|
| 267 | write(*,*) 'readwind processing wrfout file =' |
---|
| 268 | write(*,*) fnamenc |
---|
| 269 | write(*,*) 'itime, ymd, hms =', itime, jyyyymmdd, jhhmmss |
---|
| 270 | |
---|
| 271 | endif |
---|
| 272 | ! read eta_w_wrf, eta_u_wrf, p_top_wrf, ylat2d, xlon2d from the |
---|
| 273 | ! netcdf wrfout file and compare to those from the 1st met. file |
---|
| 274 | |
---|
| 275 | varname = 'ZNW' |
---|
| 276 | do i = 1, ndims_max |
---|
| 277 | lendim_exp(i) = 0 |
---|
| 278 | lendim_max(i) = 1 |
---|
| 279 | end do |
---|
| 280 | lendim_exp(1) = nwz |
---|
| 281 | lendim_max(1) = nwzmax+1 |
---|
| 282 | ndims_exp = 2 |
---|
| 283 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
| 284 | varname, dumarray_aa, & |
---|
| 285 | itime, & |
---|
| 286 | ndims, ndims_exp, ndims_max, & |
---|
| 287 | lendim, lendim_exp, lendim_max ) |
---|
| 288 | if (ierr .ne. 0) then |
---|
| 289 | fnamenc2='wrfout_d03_zn.nc' |
---|
| 290 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc2, & |
---|
| 291 | varname, dumarray_aa, & |
---|
| 292 | itime, & |
---|
| 293 | ndims, ndims_exp, ndims_max, & |
---|
| 294 | lendim, lendim_exp, lendim_max ) |
---|
| 295 | |
---|
| 296 | if (ierr .ne. 0) then |
---|
| 297 | write(*,9100) 'error doing ncread of ZNW', fnamenc |
---|
| 298 | stop |
---|
| 299 | end if |
---|
| 300 | end if |
---|
| 301 | do k = 1, nwz |
---|
| 302 | if (abs(eta_w_wrf(k) - dumarray_aa(k)) & |
---|
| 303 | .gt. toler*abs(eta_w_wrf(k))) then |
---|
| 304 | write(*,9100) 'eta_w_wrf not consistent', fnamenc |
---|
| 305 | stop |
---|
| 306 | end if |
---|
| 307 | end do |
---|
| 308 | |
---|
| 309 | varname = 'ZNU' |
---|
| 310 | lendim_exp(1) = nwz-1 |
---|
| 311 | lendim_max(1) = nwzmax+1 |
---|
| 312 | ndims_exp = 2 |
---|
| 313 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
| 314 | varname, dumarray_aa, & |
---|
| 315 | itime, & |
---|
| 316 | ndims, ndims_exp, ndims_max, & |
---|
| 317 | lendim, lendim_exp, lendim_max ) |
---|
| 318 | if (ierr .ne. 0) then |
---|
| 319 | fnamenc2='wrfout_d03_zn.nc' |
---|
| 320 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc2, & |
---|
| 321 | varname, dumarray_aa, & |
---|
| 322 | itime, & |
---|
| 323 | ndims, ndims_exp, ndims_max, & |
---|
| 324 | lendim, lendim_exp, lendim_max ) |
---|
| 325 | if (ierr .ne. 0) then |
---|
| 326 | |
---|
| 327 | write(*,9100) 'error doing ncread of ZNU', fnamenc |
---|
| 328 | stop |
---|
| 329 | end if |
---|
| 330 | end if |
---|
| 331 | do k = 1, nwz-1 |
---|
| 332 | if (abs(eta_u_wrf(k) - dumarray_aa(k)) & |
---|
| 333 | .gt. toler*abs(eta_u_wrf(k))) then |
---|
| 334 | write(*,9100) 'eta_u_wrf not consistent', fnamenc |
---|
| 335 | stop |
---|
| 336 | end if |
---|
| 337 | end do |
---|
| 338 | |
---|
| 339 | ! varname = 'P_TOP' |
---|
| 340 | ! lendim_exp(1) = 1 |
---|
| 341 | ! lendim_max(1) = 6 |
---|
| 342 | ! ndims_exp = 2 |
---|
| 343 | ! if (ext_scalar .lt. 0) ndims_exp = 1 |
---|
| 344 | ! call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
| 345 | ! varname, duma, & |
---|
| 346 | ! itime, & |
---|
| 347 | ! ndims, ndims_exp, ndims_max, & |
---|
| 348 | ! lendim, lendim_exp, lendim_max ) |
---|
| 349 | ! if (ierr .ne. 0) then |
---|
| 350 | ! write(*,9100) 'error doing ncread of P_TOP', fnamenc |
---|
| 351 | ! stop |
---|
| 352 | ! end if |
---|
| 353 | ! if (abs(p_top_wrf - duma) .gt. toler*abs(p_top_wrf)) then |
---|
| 354 | ! write(*,9100) 'p_top_wrf not consistent', fnamenc |
---|
| 355 | ! stop |
---|
| 356 | ! end if |
---|
| 357 | |
---|
| 358 | varname = 'XLAT' |
---|
| 359 | lendim_exp(1) = nx |
---|
| 360 | lendim_max(1) = nxmax |
---|
| 361 | lendim_exp(2) = ny |
---|
| 362 | lendim_max(2) = nymax |
---|
| 363 | ndims_exp = 3 |
---|
| 364 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
| 365 | varname, dumarray_pp, & |
---|
| 366 | itime, & |
---|
| 367 | ndims, ndims_exp, ndims_max, & |
---|
| 368 | lendim, lendim_exp, lendim_max ) |
---|
| 369 | if (ierr .ne. 0) then |
---|
| 370 | write(*,*) |
---|
| 371 | write(*,9100) 'error doing ncread of XLAT', fnamenc |
---|
| 372 | stop |
---|
| 373 | end if |
---|
| 374 | toler = 1.0e-6 |
---|
| 375 | do j = 0, ny-1 |
---|
| 376 | do i = 0, nx-1 |
---|
| 377 | if (abs(ylat2d(i,j) - dumarray_pp(i,j,1)) .gt. & |
---|
| 378 | toler*abs(ylat2d(i,j))) then |
---|
| 379 | write(*,9100) 'ylat2d not consistent', fnamenc |
---|
| 380 | write(*,'(a,2i5,2f16.6)') 'i,j,ylats =', i, j, & |
---|
| 381 | ylat2d(i,j), dumarray_pp(i,j,1) |
---|
| 382 | stop |
---|
| 383 | end if |
---|
| 384 | end do |
---|
| 385 | end do |
---|
| 386 | |
---|
| 387 | varname = 'XLONG' |
---|
| 388 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
| 389 | varname, dumarray_pp, & |
---|
| 390 | itime, & |
---|
| 391 | ndims, ndims_exp, ndims_max, & |
---|
| 392 | lendim, lendim_exp, lendim_max ) |
---|
| 393 | if (ierr .ne. 0) then |
---|
| 394 | write(*,*) |
---|
| 395 | write(*,9100) 'error doing ncread of XLONG', fnamenc |
---|
| 396 | stop |
---|
| 397 | end if |
---|
| 398 | do j = 0, ny-1 |
---|
| 399 | do i = 0, nx-1 |
---|
| 400 | if (abs(xlon2d(i,j) - dumarray_pp(i,j,1)) .gt. & |
---|
| 401 | toler*abs(xlon2d(i,j))) then |
---|
| 402 | write(*,9100) 'xlon2d not consistent', fnamenc |
---|
| 403 | write(*,'(a,2i5,2f16.6)') 'i,j,xlons =', i, j, & |
---|
| 404 | xlon2d(i,j), dumarray_pp(i,j,1) |
---|
| 405 | stop |
---|
| 406 | end if |
---|
| 407 | end do |
---|
| 408 | end do |
---|
| 409 | |
---|
| 410 | |
---|
| 411 | ! |
---|
| 412 | ! |
---|
| 413 | ! now read the data fields for current time |
---|
| 414 | ! the following are read from ecmwf met files |
---|
| 415 | ! U VELOCITY |
---|
| 416 | ! V VELOCITY |
---|
| 417 | ! W VELOCITY |
---|
| 418 | ! TEMPERATURE |
---|
| 419 | ! SPEC. HUMIDITY |
---|
| 420 | ! SURF. PRESS. |
---|
| 421 | ! SEA LEVEL PRESS. |
---|
| 422 | ! 10 M U VELOCITY |
---|
| 423 | ! 10 M V VELOCITY |
---|
| 424 | ! 2 M TEMPERATURE |
---|
| 425 | ! 2 M DEW POINT |
---|
| 426 | ! SNOW DEPTH |
---|
| 427 | ! CLOUD COVER |
---|
| 428 | ! LARGE SCALE PREC. |
---|
| 429 | ! CONVECTIVE PREC. |
---|
| 430 | ! SENS. HEAT FLUX |
---|
| 431 | ! SOLAR RADIATION |
---|
| 432 | ! EW SURFACE STRESS |
---|
| 433 | ! NS SURFACE STRESS |
---|
| 434 | ! ECMWF OROGRAPHY |
---|
| 435 | ! STANDARD DEVIATION OF OROGRAPHY |
---|
| 436 | ! ECMWF LAND SEA MASK |
---|
| 437 | ! |
---|
| 438 | ! |
---|
| 439 | hflswitch=.false. |
---|
| 440 | strswitch=.false. |
---|
| 441 | levdiff2=nlev_ec-nwz+1 |
---|
| 442 | |
---|
| 443 | kbgn = 1 + add_sfc_level |
---|
| 444 | ! at this point |
---|
| 445 | ! if add_sfc_level=1, then nuvz=nwz and kbgn=2 |
---|
| 446 | ! if add_sfc_level=0, then nuvz=nwz-1 and kbgn=1 |
---|
| 447 | |
---|
| 448 | ! u wind velocity |
---|
| 449 | ! the wrf output file contains (nuvz-add_sfc_level) levels |
---|
| 450 | ! read the data into k=kbgn,nuvz |
---|
| 451 | ! (interpolate it from "U-grid" to "T-grid" later) |
---|
| 452 | if (wind_option.le.0) varname = 'U' |
---|
| 453 | if (wind_option.eq.1) varname = 'AVGFLX_RUM' |
---|
| 454 | if (wind_option.eq.2) varname = 'U' |
---|
| 455 | |
---|
| 456 | do i = 1, ndims_max |
---|
| 457 | lendim_exp(i) = 0 |
---|
| 458 | lendim_max(i) = 1 |
---|
| 459 | end do |
---|
| 460 | lendim_exp(1) = nx+1 |
---|
| 461 | lendim_max(1) = nxmax |
---|
| 462 | lendim_exp(2) = ny |
---|
| 463 | lendim_max(2) = nymax |
---|
| 464 | lendim_exp(3) = nuvz-add_sfc_level |
---|
| 465 | lendim_max(3) = nuvzmax |
---|
| 466 | ndims_exp = 4 |
---|
| 467 | if (time_option.eq.0) then |
---|
| 468 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
| 469 | varname, uuh(0,0,kbgn), & |
---|
| 470 | itime, & |
---|
| 471 | ndims, ndims_exp, ndims_max, & |
---|
| 472 | lendim, lendim_exp, lendim_max ) |
---|
| 473 | if (ierr .ne. 0) then |
---|
| 474 | write(*,9100) 'error doing ncread of U', fnamenc |
---|
| 475 | if (wind_option.le.0) print*,'you asked snapshot winds' |
---|
| 476 | if (wind_option.eq.1) print*,'you asked mean winds' |
---|
| 477 | print*,'change wind_option' |
---|
| 478 | |
---|
| 479 | stop |
---|
| 480 | end if |
---|
| 481 | |
---|
| 482 | endif !test on time_option |
---|
| 483 | |
---|
| 484 | ! v wind velocity |
---|
| 485 | ! the wrf output file contains (nuvz-add_sfc_level) levels |
---|
| 486 | ! read the data into k=kbgn,nuvz |
---|
| 487 | ! (interpolate it from "V-grid" to "T-grid" later) |
---|
| 488 | if (wind_option.le.0) varname = 'V' |
---|
| 489 | if (wind_option.eq.1) varname = 'AVGFLX_RVM' |
---|
| 490 | if (wind_option.eq.2) varname = 'V' |
---|
| 491 | lendim_exp(1) = nx |
---|
| 492 | lendim_max(1) = nxmax |
---|
| 493 | lendim_exp(2) = ny+1 |
---|
| 494 | lendim_max(2) = nymax |
---|
| 495 | if (time_option.eq.0) then |
---|
| 496 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
| 497 | varname, vvh(0,0,kbgn), & |
---|
| 498 | itime, & |
---|
| 499 | ndims, ndims_exp, ndims_max, & |
---|
| 500 | lendim, lendim_exp, lendim_max ) |
---|
| 501 | if (ierr .ne. 0) then |
---|
| 502 | write(*,9100) 'error doing ncread of V', fnamenc |
---|
| 503 | if (wind_option.eq.0) print*,'you asked snapshot winds' |
---|
| 504 | if (wind_option.eq.1) print*,'you asked mean winds' |
---|
| 505 | print*,'change wind_option' |
---|
| 506 | stop |
---|
| 507 | end if |
---|
| 508 | |
---|
| 509 | endif !test on time_option |
---|
| 510 | |
---|
| 511 | ! w wind velocity |
---|
| 512 | ! this is on the "W-grid", and |
---|
| 513 | ! the wrf output file contains nwz levels, so no shifting needed |
---|
| 514 | if (wind_option.le.0) varname = 'W' |
---|
| 515 | if (wind_option.eq.1) varname = 'AVGFLX_WWM' |
---|
| 516 | if (wind_option.eq.2) varname = 'WW' |
---|
| 517 | ! print*,'varname',varname |
---|
| 518 | lendim_exp(1) = nx |
---|
| 519 | lendim_max(1) = nxmax |
---|
| 520 | lendim_exp(2) = ny |
---|
| 521 | lendim_max(2) = nymax |
---|
| 522 | lendim_exp(3) = nwz |
---|
| 523 | lendim_max(3) = nwzmax |
---|
| 524 | if (time_option.eq.0) then |
---|
| 525 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
| 526 | varname, wwh, & |
---|
| 527 | itime, & |
---|
| 528 | ndims, ndims_exp, ndims_max, & |
---|
| 529 | lendim, lendim_exp, lendim_max ) |
---|
| 530 | if (ierr .ne. 0) then |
---|
| 531 | write(*,9100) 'error doing ncread of W', fnamenc |
---|
| 532 | if (wind_option.eq.0) print*,'you asked snapshot winds' |
---|
| 533 | if (wind_option.eq.1) print*,'you asked mean winds' |
---|
| 534 | print*,'change wind_option' |
---|
| 535 | stop |
---|
| 536 | end if |
---|
| 537 | |
---|
| 538 | endif !test on time_option |
---|
| 539 | |
---|
| 540 | ! pressure - read base state and perturbation pressure, |
---|
| 541 | ! then combine |
---|
| 542 | ! the wrf output file contains (nuvz-add_sfc_level) levels |
---|
| 543 | ! read the data into k=kbgn,nuvz |
---|
| 544 | varname = 'PB' |
---|
| 545 | lendim_exp(3) = nuvz-add_sfc_level |
---|
| 546 | lendim_max(3) = nuvzmax |
---|
| 547 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
| 548 | varname, pph(0,0,kbgn,n), & |
---|
| 549 | itime, & |
---|
| 550 | ndims, ndims_exp, ndims_max, & |
---|
| 551 | lendim, lendim_exp, lendim_max ) |
---|
| 552 | if (ierr .ne. 0) then |
---|
| 553 | write(*,9100) 'error doing ncread of PB', fnamenc |
---|
| 554 | stop |
---|
| 555 | end if |
---|
| 556 | |
---|
| 557 | varname = 'P' |
---|
| 558 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
| 559 | varname, dumarray_pp(0,0,kbgn), & |
---|
| 560 | itime, & |
---|
| 561 | ndims, ndims_exp, ndims_max, & |
---|
| 562 | lendim, lendim_exp, lendim_max ) |
---|
| 563 | if (ierr .ne. 0) then |
---|
| 564 | write(*,9100) 'error doing ncread of P', fnamenc |
---|
| 565 | stop |
---|
| 566 | end if |
---|
| 567 | |
---|
| 568 | do k = kbgn, nuvz |
---|
| 569 | do j = 0, nymin1 |
---|
| 570 | do i = 0, nxmin1 |
---|
| 571 | pph(i,j,k,n) = pph(i,j,k,n) + dumarray_pp(i,j,k) |
---|
| 572 | end do |
---|
| 573 | end do |
---|
| 574 | end do |
---|
| 575 | |
---|
| 576 | |
---|
| 577 | ! height - read base state and perturbation geopotential, |
---|
| 578 | ! then combine and divide by gravity |
---|
| 579 | ! these are on the "W-grid", and |
---|
| 580 | ! the wrf output file contains nwz levels |
---|
| 581 | ! shift them also so they will be consistent with pph |
---|
| 582 | varname = 'PHB' |
---|
| 583 | lendim_exp(3) = nwz |
---|
| 584 | lendim_max(3) = nwzmax+1 |
---|
| 585 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
| 586 | varname, zzh(0,0,kbgn,n), & |
---|
| 587 | itime, & |
---|
| 588 | ndims, ndims_exp, ndims_max, & |
---|
| 589 | lendim, lendim_exp, lendim_max ) |
---|
| 590 | if (ierr .ne. 0) then |
---|
| 591 | write(*,9100) 'error doing ncread of PB', fnamenc |
---|
| 592 | stop |
---|
| 593 | end if |
---|
| 594 | |
---|
| 595 | varname = 'PH' |
---|
| 596 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
| 597 | varname, dumarray_pp(0,0,kbgn), & |
---|
| 598 | itime, & |
---|
| 599 | ndims, ndims_exp, ndims_max, & |
---|
| 600 | lendim, lendim_exp, lendim_max ) |
---|
| 601 | if (ierr .ne. 0) then |
---|
| 602 | write(*,9100) 'error doing ncread of P', fnamenc |
---|
| 603 | stop |
---|
| 604 | end if |
---|
| 605 | |
---|
| 606 | do k = kbgn, nwz+add_sfc_level |
---|
| 607 | do j = 0, nymin1 |
---|
| 608 | do i = 0, nxmin1 |
---|
| 609 | zzh(i,j,k,n) = & |
---|
| 610 | (zzh(i,j,k,n) + dumarray_pp(i,j,k))/9.81 |
---|
| 611 | end do |
---|
| 612 | end do |
---|
| 613 | end do |
---|
| 614 | |
---|
| 615 | ! now use dumarray_pp to store 1/density for stress calculation below |
---|
| 616 | ! if(sfc_option .eq. sfc_option_wrf) then |
---|
| 617 | |
---|
| 618 | ! varname = 'ALT' |
---|
| 619 | ! lendim_exp(3) = nuvz-add_sfc_level |
---|
| 620 | ! lendim_max(3) = nwzmax |
---|
| 621 | ! call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
| 622 | ! varname, dumarray_pp(0,0,kbgn), & |
---|
| 623 | ! itime, & |
---|
| 624 | ! ndims, ndims_exp, ndims_max, & |
---|
| 625 | ! lendim, lendim_exp, lendim_max ) |
---|
| 626 | ! if (ierr .ne. 0) then |
---|
| 627 | ! write(*,9100) 'error doing ncread of ALT', fnamenc |
---|
| 628 | ! stop |
---|
| 629 | ! end if |
---|
| 630 | |
---|
| 631 | ! end if |
---|
| 632 | |
---|
| 633 | ! temperature - read perturbation potential temperature, |
---|
| 634 | ! add 300. (base value), then add and convert |
---|
| 635 | ! the wrf output file contains (nuvz-add_sfc_level) levels |
---|
| 636 | ! read the data into k=kbgn,nuvz |
---|
| 637 | varname = 'T' |
---|
| 638 | lendim_exp(3) = nuvz-add_sfc_level |
---|
| 639 | lendim_max(3) = nuvzmax |
---|
| 640 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
| 641 | varname, tth(0,0,kbgn,n), & |
---|
| 642 | itime, & |
---|
| 643 | ndims, ndims_exp, ndims_max, & |
---|
| 644 | lendim, lendim_exp, lendim_max ) |
---|
| 645 | if (ierr .ne. 0) then |
---|
| 646 | write(*,9100) 'error doing ncread of T', fnamenc |
---|
| 647 | stop |
---|
| 648 | end if |
---|
| 649 | |
---|
| 650 | do k = kbgn, nuvz |
---|
| 651 | do j = 0, nymin1 |
---|
| 652 | do i = 0, nxmin1 |
---|
| 653 | ! save potential temperature to ptth |
---|
| 654 | ptth(i,j,k,n)=tth(i,j,k,n)+300. |
---|
| 655 | tth(i,j,k,n) = (tth(i,j,k,n) + 300.) * & |
---|
| 656 | (pph(i,j,k,n)/1.0e5)**0.286 |
---|
| 657 | |
---|
| 658 | end do |
---|
| 659 | end do |
---|
| 660 | end do |
---|
| 661 | |
---|
| 662 | !- |
---|
| 663 | if (turb_option .eq. turb_option_tke .or. & |
---|
| 664 | turb_option .eq. turb_option_mytke ) then |
---|
| 665 | print*, 'READ TKE',turb_option |
---|
| 666 | ! TKE - read Turbulent Kinetic, |
---|
| 667 | ! the wrf output file contains (nuvz-add_sfc_level) levels |
---|
| 668 | ! read the data into k=kbgn,nuvz |
---|
| 669 | varname = 'TKE' |
---|
| 670 | lendim_exp(3) = nuvz-add_sfc_level |
---|
| 671 | lendim_max(3) = nuvzmax |
---|
| 672 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
| 673 | varname, tkeh(0,0,kbgn,n), & |
---|
| 674 | itime, & |
---|
| 675 | ndims, ndims_exp, ndims_max, & |
---|
| 676 | lendim, lendim_exp, lendim_max ) |
---|
| 677 | |
---|
| 678 | if (ierr .ne. 0) then |
---|
| 679 | print*,'NO TKE available. Try TKE_PBL instead' |
---|
| 680 | varname = 'TKE_PBL' |
---|
| 681 | lendim_exp(3) = nuvz-add_sfc_level |
---|
| 682 | lendim_max(3) = nuvzmax |
---|
| 683 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
| 684 | varname, tkeh(0,0,kbgn,n), & |
---|
| 685 | itime, & |
---|
| 686 | ndims, ndims_exp, ndims_max, & |
---|
| 687 | lendim, lendim_exp, lendim_max ) |
---|
| 688 | endif |
---|
| 689 | if (ierr .ne. 0) then |
---|
| 690 | print*,'NO TKE_PBL available. Try QKE instead' |
---|
| 691 | varname = 'qke' |
---|
| 692 | lendim_exp(3) = nuvz-add_sfc_level |
---|
| 693 | lendim_max(3) = nuvzmax |
---|
| 694 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
| 695 | varname, tkeh(0,0,kbgn,n), & |
---|
| 696 | itime, & |
---|
| 697 | ndims, ndims_exp, ndims_max, & |
---|
| 698 | lendim, lendim_exp, lendim_max ) |
---|
| 699 | tkeh=tkeh/2. !conversion of qke |
---|
| 700 | endif |
---|
| 701 | if (ierr .ne. 0) then |
---|
| 702 | varname = 'QKE' |
---|
| 703 | lendim_exp(3) = nuvz-add_sfc_level |
---|
| 704 | lendim_max(3) = nuvzmax |
---|
| 705 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
| 706 | varname, tkeh(0,0,kbgn,n), & |
---|
| 707 | itime, & |
---|
| 708 | ndims, ndims_exp, ndims_max, & |
---|
| 709 | lendim, lendim_exp, lendim_max ) |
---|
| 710 | tkeh=tkeh/2. !conversion of qke |
---|
| 711 | endif |
---|
| 712 | |
---|
| 713 | if (ierr .ne. 0) then |
---|
| 714 | write(*,9100) 'error doing ncread of TKE', fnamenc |
---|
| 715 | write(*,*)'change turb_option NOT to use TKE or change input file' |
---|
| 716 | print*,'change SFC_OPTION to 0 as well' |
---|
| 717 | stop |
---|
| 718 | end if |
---|
| 719 | |
---|
| 720 | endif |
---|
| 721 | |
---|
| 722 | |
---|
| 723 | !- |
---|
| 724 | |
---|
| 725 | ! specific humidity - read mixing ratio (kg-water-vapor/kg-dry-air), |
---|
| 726 | ! then convert to (kg-water-vapor/kg-moist-air) |
---|
| 727 | ! the wrf output file contains (nuvz-add_sfc_level) levels |
---|
| 728 | ! read the data into k=kbgn,nuvz |
---|
| 729 | varname = 'QVAPOR' |
---|
| 730 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
| 731 | varname, qvh(0,0,kbgn,n), & |
---|
| 732 | itime, & |
---|
| 733 | ndims, ndims_exp, ndims_max, & |
---|
| 734 | lendim, lendim_exp, lendim_max ) |
---|
| 735 | if (ierr .ne. 0) then |
---|
| 736 | write(*,9100) 'error doing ncread of QVAPOR', fnamenc |
---|
| 737 | stop |
---|
| 738 | end if |
---|
| 739 | |
---|
| 740 | do k = kbgn, nuvz |
---|
| 741 | do j = 0, nymin1 |
---|
| 742 | do i = 0, nxmin1 |
---|
| 743 | qvh(i,j,k,n) = max( qvh(i,j,k,n), 0.0 ) |
---|
| 744 | qvh(i,j,k,n) = qvh(i,j,k,n)/(1.0 + qvh(i,j,k,n)) |
---|
| 745 | end do |
---|
| 746 | end do |
---|
| 747 | end do |
---|
| 748 | |
---|
| 749 | |
---|
| 750 | ! surface pressure |
---|
| 751 | varname = 'PSFC' |
---|
| 752 | lendim_exp(3) = 0 |
---|
| 753 | lendim_max(3) = 1 |
---|
| 754 | ndims_exp = 3 |
---|
| 755 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
| 756 | varname, ps(0,0,1,n), & |
---|
| 757 | itime, & |
---|
| 758 | ndims, ndims_exp, ndims_max, & |
---|
| 759 | lendim, lendim_exp, lendim_max ) |
---|
| 760 | if (ierr .ne. 0) then |
---|
| 761 | write(*,9100) 'error doing ncread of PSFC', fnamenc |
---|
| 762 | stop |
---|
| 763 | end if |
---|
| 764 | |
---|
| 765 | ! for the mexico city grid 3 simulation, the surface and |
---|
| 766 | ! level 1 pressures are not as consistent as one would like, |
---|
| 767 | ! with the problems occuring near the domain boundaries. |
---|
| 768 | ! so do the following |
---|
| 769 | ! -- calculate surface pressure from lowest level pressure, temp, height |
---|
| 770 | ! -- use wrf pressures (pph array) wherever possible |
---|
| 771 | ! (avoid using surface pressure and the akz/bkz, akm/bkm) |
---|
| 772 | do j = 0, nymin1 |
---|
| 773 | do i = 0, nxmin1 |
---|
| 774 | duma = ps(i,j,1,n) |
---|
| 775 | dumdz = 0.5*(zzh(i,j,kbgn+1,n) - zzh(i,j,kbgn,n)) |
---|
| 776 | tv = tth(i,j,kbgn,n)*(1.+0.61*qvh(i,j,kbgn,n)) |
---|
| 777 | ps(i,j,1,n) = pph(i,j,kbgn,n)*exp( dumdz*ga/(r_air*tv) ) |
---|
| 778 | end do |
---|
| 779 | end do |
---|
| 780 | |
---|
| 781 | |
---|
| 782 | ! 10 meter u velocity |
---|
| 783 | ! note: u10 is on the "T-grid" already |
---|
| 784 | varname = 'U10' |
---|
| 785 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
| 786 | varname, u10(0,0,1,n), & |
---|
| 787 | itime, & |
---|
| 788 | ndims, ndims_exp, ndims_max, & |
---|
| 789 | lendim, lendim_exp, lendim_max ) |
---|
| 790 | if (ierr .ne. 0) then |
---|
| 791 | write(*,9100) 'error doing ncread of U10', fnamenc |
---|
| 792 | stop |
---|
| 793 | end if |
---|
| 794 | |
---|
| 795 | |
---|
| 796 | ! 10 meter v velocity |
---|
| 797 | ! note: v10 is on the "T-grid" already |
---|
| 798 | varname = 'V10' |
---|
| 799 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
| 800 | varname, v10(0,0,1,n), & |
---|
| 801 | itime, & |
---|
| 802 | ndims, ndims_exp, ndims_max, & |
---|
| 803 | lendim, lendim_exp, lendim_max ) |
---|
| 804 | if (ierr .ne. 0) then |
---|
| 805 | write(*,9100) 'error doing ncread of V10', fnamenc |
---|
| 806 | stop |
---|
| 807 | end if |
---|
| 808 | |
---|
| 809 | |
---|
| 810 | ! 2 meter temperature |
---|
| 811 | varname = 'T2' |
---|
| 812 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
| 813 | varname, tt2(0,0,1,n), & |
---|
| 814 | itime, & |
---|
| 815 | ndims, ndims_exp, ndims_max, & |
---|
| 816 | lendim, lendim_exp, lendim_max ) |
---|
| 817 | if (ierr .ne. 0) then |
---|
| 818 | write(*,9100) 'error doing ncread of T2', fnamenc |
---|
| 819 | stop |
---|
| 820 | end if |
---|
| 821 | |
---|
| 822 | |
---|
| 823 | ! 2 meter dew point - read 2 meter water vapor mixing ratio |
---|
| 824 | ! then calculate the dew point |
---|
| 825 | varname = 'Q2' |
---|
| 826 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
| 827 | varname, td2(0,0,1,n), & |
---|
| 828 | itime, & |
---|
| 829 | ndims, ndims_exp, ndims_max, & |
---|
| 830 | lendim, lendim_exp, lendim_max ) |
---|
| 831 | if (ierr .ne. 0) then |
---|
| 832 | write(*,9100) 'error doing ncread of Q2' |
---|
| 833 | !, fnamenc |
---|
| 834 | do j = 0, nymin1 |
---|
| 835 | do i = 0, nxmin1 |
---|
| 836 | ! 29-nov-2005 - changed qvh(i,j,1,n) to qvh(i,j,kbgn,n) here |
---|
| 837 | td2(i,j,1,n) = qvh(i,j,kbgn,n) |
---|
| 838 | end do |
---|
| 839 | end do |
---|
| 840 | end if |
---|
| 841 | |
---|
| 842 | if (wind_option.ge.1) then |
---|
| 843 | ! print*,'mean wind from WRF is used' |
---|
| 844 | ! print*,'option ',wind_option |
---|
| 845 | varname = 'MU' |
---|
| 846 | |
---|
| 847 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
| 848 | varname, mu(0,0,1), & |
---|
| 849 | itime, & |
---|
| 850 | ndims, ndims_exp, ndims_max, & |
---|
| 851 | lendim, lendim_exp, lendim_max ) |
---|
| 852 | if (ierr .ne. 0) then |
---|
| 853 | write(*,9100) 'error doing MU', fnamenc |
---|
| 854 | stop |
---|
| 855 | end if |
---|
| 856 | |
---|
| 857 | varname = 'MUB' |
---|
| 858 | |
---|
| 859 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
| 860 | varname, mub(0,0,1), & |
---|
| 861 | itime, & |
---|
| 862 | ndims, ndims_exp, ndims_max, & |
---|
| 863 | lendim, lendim_exp, lendim_max ) |
---|
| 864 | if (ierr .ne. 0) then |
---|
| 865 | write(*,9100) 'error doing MUV', fnamenc |
---|
| 866 | stop |
---|
| 867 | end if |
---|
| 868 | endif |
---|
| 869 | |
---|
| 870 | ! varname = 'MAPFAC_UX' |
---|
| 871 | ! varname = 'MAPFAC_UY' !try |
---|
| 872 | ! lendim_exp(1) = nx+1 |
---|
| 873 | ! lendim_max(1) = nxmax |
---|
| 874 | ! lendim_exp(2) = ny |
---|
| 875 | ! lendim_max(2) = nymax |
---|
| 876 | ! |
---|
| 877 | ! call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
| 878 | ! varname, m_u(0,0,1), & |
---|
| 879 | ! itime, & |
---|
| 880 | ! ndims, ndims_exp, ndims_max, & |
---|
| 881 | ! lendim, lendim_exp, lendim_max ) |
---|
| 882 | ! if (ierr .ne. 0) then |
---|
| 883 | ! write(*,9100) 'error doing MAP U', fnamenc |
---|
| 884 | ! stop |
---|
| 885 | ! end if |
---|
| 886 | ! |
---|
| 887 | ! varname = 'MAPFAC_VY' |
---|
| 888 | ! varname = 'MAPFAC_VX' !try |
---|
| 889 | ! lendim_exp(1) = nx |
---|
| 890 | ! lendim_max(1) = nxmax |
---|
| 891 | ! lendim_exp(2) = ny+1 |
---|
| 892 | ! lendim_max(2) = nymax |
---|
| 893 | ! |
---|
| 894 | ! call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
| 895 | ! varname, m_v(0,0,1), & |
---|
| 896 | ! itime, & |
---|
| 897 | ! ndims, ndims_exp, ndims_max, & |
---|
| 898 | ! lendim, lendim_exp, lendim_max ) |
---|
| 899 | ! if (ierr .ne. 0) then |
---|
| 900 | ! write(*,9100) 'error doing MAP V', fnamenc |
---|
| 901 | ! stop |
---|
| 902 | ! end if |
---|
| 903 | |
---|
| 904 | ! varname = 'MAPFAC_MX' |
---|
| 905 | ! lendim_exp(1) = nx |
---|
| 906 | ! lendim_max(1) = nxmax |
---|
| 907 | ! lendim_exp(2) = ny |
---|
| 908 | ! lendim_max(2) = nymax |
---|
| 909 | ! |
---|
| 910 | ! call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
| 911 | ! varname, m_x(0,0,1), & |
---|
| 912 | ! itime, & |
---|
| 913 | ! ndims, ndims_exp, ndims_max, & |
---|
| 914 | ! lendim, lendim_exp, lendim_max ) |
---|
| 915 | ! if (ierr .ne. 0) then |
---|
| 916 | ! write(*,9100) 'error doing MAP X', fnamenc |
---|
| 917 | ! varname = 'MAPFAC_U' |
---|
| 918 | ! lendim_exp(1) = nx+1 |
---|
| 919 | ! lendim_max(1) = nxmax |
---|
| 920 | ! lendim_exp(2) = ny |
---|
| 921 | ! lendim_max(2) = nymax |
---|
| 922 | ! call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
| 923 | ! varname, m_u(0,0,1), & |
---|
| 924 | ! itime, & |
---|
| 925 | ! ndims, ndims_exp, ndims_max, & |
---|
| 926 | ! lendim, lendim_exp, lendim_max ) |
---|
| 927 | ! do j = 0, nymin1 |
---|
| 928 | ! do i = 0, nxmin1 |
---|
| 929 | ! m_x(i,j,1)=(m_u(i,j,1)+m_u(i+1,j,1))*0.5 |
---|
| 930 | ! enddo |
---|
| 931 | ! enddo |
---|
| 932 | ! if (ierr .ne. 0) then |
---|
| 933 | ! write(*,9100) 'error doing MAP U', fnamenc |
---|
| 934 | ! print*,'NO MAP FACTOR IS GOING TO BE USED.' |
---|
| 935 | ! print*,'LARGE UNCERTAINTIES TO BE EXPECTED' |
---|
| 936 | ! do j = 0, nymin1 |
---|
| 937 | ! do i = 0, nxmin1 |
---|
| 938 | ! m_x(i,j,1)=1. |
---|
| 939 | ! enddo |
---|
| 940 | ! enddo |
---|
| 941 | ! end if |
---|
| 942 | ! end if |
---|
| 943 | |
---|
| 944 | ! varname = 'MAPFAC_M' |
---|
| 945 | ! lendim_exp(1) = nx |
---|
| 946 | ! lendim_max(1) = nxmax |
---|
| 947 | ! lendim_exp(2) = ny |
---|
| 948 | ! lendim_max(2) = nymax |
---|
| 949 | ! |
---|
| 950 | ! call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
| 951 | ! varname, m_z(0,0,1), & |
---|
| 952 | ! itime, & |
---|
| 953 | ! ndims, ndims_exp, ndims_max, & |
---|
| 954 | ! lendim, lendim_exp, lendim_max ) |
---|
| 955 | ! if (ierr .ne. 0) then |
---|
| 956 | ! write(*,9100) 'error doing MAP W', fnamenc |
---|
| 957 | ! stop |
---|
| 958 | ! end if |
---|
| 959 | |
---|
| 960 | ! varname = 'MAPFAC_MY' |
---|
| 961 | ! lendim_exp(1) = nx |
---|
| 962 | ! lendim_max(1) = nxmax |
---|
| 963 | ! lendim_exp(2) = ny |
---|
| 964 | ! lendim_max(2) = nymax |
---|
| 965 | ! |
---|
| 966 | ! call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
| 967 | ! varname, m_y(0,0,1), & |
---|
| 968 | ! itime, & |
---|
| 969 | ! ndims, ndims_exp, ndims_max, & |
---|
| 970 | ! lendim, lendim_exp, lendim_max ) |
---|
| 971 | ! if (ierr .ne. 0) then |
---|
| 972 | ! write(*,9100) 'error doing MAP Y', fnamenc |
---|
| 973 | ! varname = 'MAPFAC_V' |
---|
| 974 | ! lendim_exp(1) = nx |
---|
| 975 | ! lendim_max(1) = nxmax |
---|
| 976 | ! lendim_exp(2) = ny+1 |
---|
| 977 | ! lendim_max(2) = nymax |
---|
| 978 | ! call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
| 979 | ! varname, m_v(0,0,1), & |
---|
| 980 | ! itime, & |
---|
| 981 | ! ndims, ndims_exp, ndims_max, & |
---|
| 982 | ! lendim, lendim_exp, lendim_max ) |
---|
| 983 | ! do j = 0, nymin1 |
---|
| 984 | ! do i = 0, nxmin1 |
---|
| 985 | ! m_y(i,j,1)=(m_v(i,j,1)+m_v(i,j+1,1))*0.5 |
---|
| 986 | ! enddo |
---|
| 987 | ! enddo |
---|
| 988 | ! if (ierr .ne. 0) then |
---|
| 989 | ! write(*,9100) 'ERROR doing MAP V', fnamenc |
---|
| 990 | ! print*,'NO MAP FACTOR IS GOING TO BE USED.' |
---|
| 991 | ! print*,'LARGE UNCERTAINTIES TO BE EXPECTED' |
---|
| 992 | ! do j = 0, nymin1 |
---|
| 993 | ! do i = 0, nxmin1 |
---|
| 994 | ! m_y(i,j,1)=1. |
---|
| 995 | ! enddo |
---|
| 996 | ! enddo |
---|
| 997 | ! end if |
---|
| 998 | ! end if |
---|
| 999 | lendim_exp(1) = nx |
---|
| 1000 | lendim_max(1) = nxmax |
---|
| 1001 | lendim_exp(2) = ny |
---|
| 1002 | lendim_max(2) = nymax |
---|
| 1003 | |
---|
| 1004 | |
---|
| 1005 | |
---|
| 1006 | ! calculate water vapor pressure in mb, from sfc pressure |
---|
| 1007 | ! and 2 m mixing ratio |
---|
| 1008 | iduma = 0 |
---|
| 1009 | do j = 0, nymin1 |
---|
| 1010 | do i = 0, nxmin1 |
---|
| 1011 | ! 29-nov-2005 - added this to catch occasional tt2n=0.0 values |
---|
| 1012 | duma = max( 100.0, tth(i,j,kbgn,n)-50.0 ) |
---|
| 1013 | if (tt2(i,j,1,n) .le. duma) then |
---|
| 1014 | iduma = iduma + 1 |
---|
| 1015 | if (iduma .eq. 1) then |
---|
| 1016 | write(*,*) 'readwind - bad tt2 at' |
---|
| 1017 | write(*,*) 'i, j, tt2 =', i, j, tt2(i,j,1,n) |
---|
| 1018 | end if |
---|
| 1019 | ! stop |
---|
| 1020 | tt2(i,j,1,n) = tth(i,j,kbgn,n) |
---|
| 1021 | td2(i,j,1,n) = qvh(i,j,kbgn,n) |
---|
| 1022 | end if |
---|
| 1023 | duma = td2(i,j,1,n)/0.622 |
---|
| 1024 | ewater_mb = 0.01*( 0.99976*ps(i,j,1,n)*duma/(1.0+duma) ) |
---|
| 1025 | esatwater_mb = 0.01*ew(tt2(i,j,1,n)) |
---|
| 1026 | ewater_mb = max( 1.0e-10, min( esatwater_mb, ewater_mb ) ) |
---|
| 1027 | ! then use the following, which is from an old 1970's report |
---|
| 1028 | ! (reference not available, but the formula works) |
---|
| 1029 | ! tdew(in C) = (4318.76/(19.5166 - ln(ewater(in mb)))) - 243.893 |
---|
| 1030 | td2(i,j,1,n) = 273.16 + & |
---|
| 1031 | (4318.76/(19.5166 - log(ewater_mb))) - 243.893 |
---|
| 1032 | end do |
---|
| 1033 | end do |
---|
| 1034 | if (iduma .gt. 0) write(*,*) & |
---|
| 1035 | 'readwind - bad tt2 count =', iduma |
---|
| 1036 | |
---|
| 1037 | |
---|
| 1038 | ! sea level pressure - calculate it from surface pressure and |
---|
| 1039 | ! ground elevation using standard atmosphere relations |
---|
| 1040 | do j = 0, nymin1 |
---|
| 1041 | do i = 0, nxmin1 |
---|
| 1042 | msl(i,j,1,n) = ps(i,j,1,n)/ & |
---|
| 1043 | ((1.0 - 6.5e-3*oro(i,j)/288.0)**5.2553) |
---|
| 1044 | end do |
---|
| 1045 | end do |
---|
| 1046 | |
---|
| 1047 | |
---|
| 1048 | ! large scale precipitation |
---|
| 1049 | ! convective precipitation |
---|
| 1050 | ! the wrf output files contain these as "accumulated totals" |
---|
| 1051 | ! I need to find out if these are accumulated over the output |
---|
| 1052 | ! file frequency, or over the total run. |
---|
| 1053 | ! For now, set to zero |
---|
| 1054 | ! total cloud cover |
---|
| 1055 | ! Doesn't appear to be any 2-d cloud cover field in the |
---|
| 1056 | ! wrf output. |
---|
| 1057 | ! For now, set to zero |
---|
| 1058 | do j = 0, nymin1 |
---|
| 1059 | do i = 0, nxmin1 |
---|
| 1060 | lsprec(i,j,1,n) = 0.0 |
---|
| 1061 | convprec(i,j,1,n) = 0.0 |
---|
| 1062 | tcc(i,j,1,n) = 0.0 |
---|
| 1063 | end do |
---|
| 1064 | end do |
---|
| 1065 | |
---|
| 1066 | ! |
---|
| 1067 | ! Large-scale precipitation (accumulated rain, mm) |
---|
| 1068 | ! will convert to mm/h in interpolat_rain.f |
---|
| 1069 | varname = 'RAINNC' |
---|
| 1070 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
| 1071 | varname, lsprec(0,0,1,n), & |
---|
| 1072 | itime, & |
---|
| 1073 | dims, ndims_exp, ndims_max, & |
---|
| 1074 | lendim, lendim_exp, lendim_max ) |
---|
| 1075 | if (ierr .ne. 0) then |
---|
| 1076 | write(*,9100) 'error doing ncread of RAINNC,set to zero', fnamenc |
---|
| 1077 | do j = 0, nymin1 |
---|
| 1078 | do i = 0, nxmin1 |
---|
| 1079 | lsprec(i,j,1,n) = 0.0 |
---|
| 1080 | end do |
---|
| 1081 | end do |
---|
| 1082 | end if |
---|
| 1083 | |
---|
| 1084 | ! |
---|
| 1085 | ! Convective cumulus precipitation (accumulated rain, mm) |
---|
| 1086 | |
---|
| 1087 | varname = 'RAINC' |
---|
| 1088 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
| 1089 | varname, convprec(0,0,1,n), & |
---|
| 1090 | itime, & |
---|
| 1091 | ndims, ndims_exp, ndims_max, & |
---|
| 1092 | lendim, lendim_exp, lendim_max ) |
---|
| 1093 | if (ierr .ne. 0) then |
---|
| 1094 | write(*,9100) 'error doing ncread of RAINC, set to zero', fnamenc |
---|
| 1095 | do j = 0, nymin1 |
---|
| 1096 | do i = 0, nxmin1 |
---|
| 1097 | convprec(i,j,1,n) = 0.0 |
---|
| 1098 | end do |
---|
| 1099 | end do |
---|
| 1100 | end if |
---|
| 1101 | |
---|
| 1102 | ! |
---|
| 1103 | ! Clound fraction (cloud cover) |
---|
| 1104 | ! varname = 'CLDFRA' |
---|
| 1105 | ! call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
| 1106 | ! varname, tcc(0,0,1,n), & |
---|
| 1107 | ! itime, & |
---|
| 1108 | ! ndims, ndims_exp, ndims_max, & |
---|
| 1109 | ! lendim, lendim_exp, lendim_max ) |
---|
| 1110 | ! if (ierr .ne. 0) then |
---|
| 1111 | !! write(*,9100) 'error doing ncread of CLDFRA, set to zero' |
---|
| 1112 | !!, fnamenc |
---|
| 1113 | ! do j = 0, nymin1 |
---|
| 1114 | ! do i = 0, nxmin1 |
---|
| 1115 | ! tcc(i,j,1,n) = 0.0 |
---|
| 1116 | ! end do |
---|
| 1117 | ! end do |
---|
| 1118 | ! end if |
---|
| 1119 | !!C write(*,*)'read CLDFRA 0-sucess ',ierr |
---|
| 1120 | |
---|
| 1121 | ! |
---|
| 1122 | |
---|
| 1123 | ! snow depth |
---|
| 1124 | varname = 'SNOWH' |
---|
| 1125 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
| 1126 | varname, sd(0,0,1,n), & |
---|
| 1127 | itime, & |
---|
| 1128 | ndims, ndims_exp, ndims_max, & |
---|
| 1129 | lendim, lendim_exp, lendim_max ) |
---|
| 1130 | if (ierr .ne. 0) then |
---|
| 1131 | ! write(*,9100) 'error doing ncread of SNOWH', fnamenc |
---|
| 1132 | do j = 0, nymin1 |
---|
| 1133 | do i = 0, nxmin1 |
---|
| 1134 | sd(i,j,1,n) = 0.0 |
---|
| 1135 | end do |
---|
| 1136 | end do |
---|
| 1137 | end if |
---|
| 1138 | |
---|
| 1139 | |
---|
| 1140 | ! surface sensible heat flux (positive <--> upwards) |
---|
| 1141 | varname = 'HFX' |
---|
| 1142 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
| 1143 | varname, sshf(0,0,1,n), & |
---|
| 1144 | itime, & |
---|
| 1145 | ndims, ndims_exp, ndims_max, & |
---|
| 1146 | lendim, lendim_exp, lendim_max ) |
---|
| 1147 | do j = 0, nymin1 |
---|
| 1148 | do i = 0, nxmin1 |
---|
| 1149 | sshf(i,j,1,n) = -sshf(i,j,1,n) |
---|
| 1150 | end do |
---|
| 1151 | end do |
---|
| 1152 | |
---|
| 1153 | if (ierr .ne. 0) then |
---|
| 1154 | write(*,9100) 'error doing ncread of HFX', fnamenc |
---|
| 1155 | do j = 0, nymin1 |
---|
| 1156 | do i = 0, nxmin1 |
---|
| 1157 | sshf(i,j,1,n) = 0.0 |
---|
| 1158 | end do |
---|
| 1159 | end do |
---|
| 1160 | hflswitch=.false. ! Heat flux is not available |
---|
| 1161 | else |
---|
| 1162 | hflswitch=.true. ! Heat flux is available |
---|
| 1163 | ! limit to values to bounds originally used by flexpart? |
---|
| 1164 | ! do 1502 j=0,nymin1 |
---|
| 1165 | ! do 1502 i=0,nxmin1 |
---|
| 1166 | ! if(sshf(i,j,1,n).gt.200.) sshf(i,j,1,n)=200. |
---|
| 1167 | ! if(sshf(i,j,1,n).lt.-400.) sshf(i,j,1,n)=-400. |
---|
| 1168 | !1502 continue |
---|
| 1169 | end if |
---|
| 1170 | |
---|
| 1171 | ! ustar |
---|
| 1172 | varname = 'UST' |
---|
| 1173 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
| 1174 | varname, ustar(0,0,1,n), & |
---|
| 1175 | itime, & |
---|
| 1176 | ndims, ndims_exp, ndims_max, & |
---|
| 1177 | lendim, lendim_exp, lendim_max ) |
---|
| 1178 | if (ierr .ne. 0) then |
---|
| 1179 | write(*,9100) 'error doing ncread of UST', fnamenc |
---|
| 1180 | do j = 0, nymin1 |
---|
| 1181 | do i = 0, nxmin1 |
---|
| 1182 | ustar(i,j,1,n) = 0.0 |
---|
| 1183 | end do |
---|
| 1184 | end do |
---|
| 1185 | strswitch=.false. ! ustar is not available |
---|
| 1186 | else |
---|
| 1187 | strswitch=.true. ! ustar is available |
---|
| 1188 | do j=0,nymin1 |
---|
| 1189 | do i=0,nxmin1 |
---|
| 1190 | ! surfstr(i,j,1,n)=ustar(i,j,1,n)/dumarray_pp(i,j,kbgn) |
---|
| 1191 | enddo |
---|
| 1192 | enddo |
---|
| 1193 | end if |
---|
| 1194 | |
---|
| 1195 | ! pblh |
---|
| 1196 | if(sfc_option .eq. sfc_option_wrf) then |
---|
| 1197 | varname = 'PBLH' |
---|
| 1198 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
| 1199 | varname, hmix(0,0,1,n), & |
---|
| 1200 | itime, & |
---|
| 1201 | ndims, ndims_exp, ndims_max, & |
---|
| 1202 | lendim, lendim_exp, lendim_max ) |
---|
| 1203 | if (ierr .ne. 0) then |
---|
| 1204 | write(*,9100) 'error doing ncread of PBLH', fnamenc |
---|
| 1205 | stop |
---|
| 1206 | endif |
---|
| 1207 | |
---|
| 1208 | endif |
---|
| 1209 | |
---|
| 1210 | ! surface solar radiation flux (positive <--> downwards) |
---|
| 1211 | varname = 'SWDOWN' |
---|
| 1212 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
| 1213 | varname, ssr(0,0,1,n), & |
---|
| 1214 | itime, & |
---|
| 1215 | ndims, ndims_exp, ndims_max, & |
---|
| 1216 | lendim, lendim_exp, lendim_max ) |
---|
| 1217 | if (ierr .ne. 0) then |
---|
| 1218 | write(*,9100) 'error doing ncread of SWDOWN', fnamenc |
---|
| 1219 | do j = 0, nymin1 |
---|
| 1220 | do i = 0, nxmin1 |
---|
| 1221 | ssr(i,j,1,n) = 0.0 |
---|
| 1222 | end do |
---|
| 1223 | end do |
---|
| 1224 | else |
---|
| 1225 | do j = 0, nymin1 |
---|
| 1226 | do i = 0, nxmin1 |
---|
| 1227 | ssr(i,j,1,n) = max( ssr(i,j,1,n), 0.0 ) |
---|
| 1228 | end do |
---|
| 1229 | end do |
---|
| 1230 | end if |
---|
| 1231 | |
---|
| 1232 | |
---|
| 1233 | ! ew & ns surface stress |
---|
| 1234 | ! Doesn't appear to be any 2-d cloud cover field in the |
---|
| 1235 | ! wrf output. |
---|
| 1236 | ! For now, set to zero |
---|
| 1237 | do j = 0, nymin1 |
---|
| 1238 | do i = 0, nxmin1 |
---|
| 1239 | ewss(i,j) = 0.0 |
---|
| 1240 | nsss(i,j) = 0.0 |
---|
| 1241 | end do |
---|
| 1242 | end do |
---|
| 1243 | ! strswitch=.false. ! Surface stress is not available |
---|
| 1244 | |
---|
| 1245 | |
---|
| 1246 | ! orography |
---|
| 1247 | ! standard deviation of orography |
---|
| 1248 | ! land sea mask |
---|
| 1249 | ! these should be fixed during a simulation |
---|
| 1250 | ! so there is no reason to do them again ?? |
---|
| 1251 | |
---|
| 1252 | |
---|
| 1253 | ! *** done with reading the wrf output file *** |
---|
| 1254 | |
---|
| 1255 | ! print*,'uu out1',uuh(0,259,1:10) |
---|
| 1256 | ! print*,'mu out1',mu(0,259,1),mub(0,259,1) |
---|
| 1257 | ! print*,'m_xn out1',m_x(0,259,1),m_y(0,259,1) |
---|
| 1258 | |
---|
| 1259 | |
---|
| 1260 | ! interpolate uuh from the "U-grid" to the "T-grid" |
---|
| 1261 | ! interpolate vvh from the "V-grid" to the "T-grid" |
---|
| 1262 | ! new: convert mass weighted wind to wind. |
---|
| 1263 | ! print*,'wind_option',wind_option |
---|
| 1264 | if (wind_option.le.0) then |
---|
| 1265 | ! if (wind_option.eq.-1) then |
---|
| 1266 | ! call calc_uvmet(uuh,vvh,urot,vrot,1) |
---|
| 1267 | ! endif |
---|
| 1268 | do k = kbgn, nuvz |
---|
| 1269 | do j = 0, nymin1 |
---|
| 1270 | do i = 0, nxmin1 |
---|
| 1271 | ! needs to rotate u and v to work. needs read alpha or something like that. |
---|
| 1272 | ! if in mercator, no need. |
---|
| 1273 | ! divh(i,j,k)=(uuh(i+1,j,k)-uuh(i,j,k))/dx & |
---|
| 1274 | ! +(vvh(i,j+1,k)-vvh(i,j,k))/dy |
---|
| 1275 | if (wind_option.lt.0) then |
---|
| 1276 | divh(i,j,k)=(uuh(i+1,j,k)-uuh(i,j,k))/dx*m_x(i,j,1) & |
---|
| 1277 | +(vvh(i,j+1,k)-vvh(i,j,k))/dy*m_y(i,j,1) |
---|
| 1278 | endif |
---|
| 1279 | uuh(i,j,k) = 0.5*(uuh(i,j,k) + uuh(i+1,j,k)) |
---|
| 1280 | vvh(i,j,k) = 0.5*(vvh(i,j,k) + vvh(i,j+1,k)) |
---|
| 1281 | end do |
---|
| 1282 | end do |
---|
| 1283 | end do |
---|
| 1284 | elseif (wind_option.eq.1) then |
---|
| 1285 | do k = kbgn, nuvz |
---|
| 1286 | do j = 0, nymin1 |
---|
| 1287 | do i = 0, nxmin1 |
---|
| 1288 | ! uuh(i,j,k) = 0.5*(uuh(i,j,k)*m_u(i,j,1) + uuh(i+1,j,k)*m_u(i+1,j,1)) |
---|
| 1289 | ! vvh(i,j,k) = 0.5*(vvh(i,j,k)*m_v(i,j,1) + vvh(i,j+1,k)*m_v(i,j+1,1)) |
---|
| 1290 | uuh(i,j,k) = 0.5*(uuh(i,j,k) + uuh(i+1,j,k)) |
---|
| 1291 | vvh(i,j,k) = 0.5*(vvh(i,j,k) + vvh(i,j+1,k)) |
---|
| 1292 | mu2=mu(i,j,1)+mub(i,j,1) |
---|
| 1293 | ! uuh(i,j,k) = uuh(i,j,k)/mu2!*0.5*(m_u(i,j,1)+m_u(i+1,j,1)) |
---|
| 1294 | ! vvh(i,j,k) = vvh(i,j,k)/mu2!*0.5*(m_v(i,j,1)+m_v(i,j+1,1)) |
---|
| 1295 | ! m_y=0.5*(m_v(i,j,1)+m_v(i,j+1,1)) |
---|
| 1296 | uuh(i,j,k) = uuh(i,j,k)/mu2 !*m_y(i,j,1) !without m=true wind |
---|
| 1297 | vvh(i,j,k) = vvh(i,j,k)/mu2 !*m_x(i,j,1) |
---|
| 1298 | wwh(i,j,k) = wwh(i,j,k)/mu2 !*m_y(i,j,1) |
---|
| 1299 | end do |
---|
| 1300 | end do |
---|
| 1301 | end do |
---|
| 1302 | elseif (wind_option.eq.2) then |
---|
| 1303 | do k = kbgn, nuvz |
---|
| 1304 | do j = 0, nymin1 |
---|
| 1305 | do i = 0, nxmin1 |
---|
| 1306 | uuh(i,j,k) = 0.5*(uuh(i,j,k) + uuh(i+1,j,k)) |
---|
| 1307 | vvh(i,j,k) = 0.5*(vvh(i,j,k) + vvh(i,j+1,k)) |
---|
| 1308 | mu2=mu(i,j,1)+mub(i,j,1) |
---|
| 1309 | wwh(i,j,k) = wwh(i,j,k)/mu2 !*m_y(i,j,1) |
---|
| 1310 | end do |
---|
| 1311 | end do |
---|
| 1312 | end do |
---|
| 1313 | |
---|
| 1314 | endif |
---|
| 1315 | |
---|
| 1316 | |
---|
| 1317 | ! for ecmwf flexpart, if nwz = nlev_ec+1, then wwh is set |
---|
| 1318 | ! to zero at the top level |
---|
| 1319 | ! for wrf, nlev_ec==n_bottom_top and so nwz = nlev_ec+1. |
---|
| 1320 | ! however, it doesn't seem appropriate to zero wwh at |
---|
| 1321 | ! the model top which might be ~100 hPa. |
---|
| 1322 | ! so deactivate this for now |
---|
| 1323 | ! if(levdiff2.eq.0) then |
---|
| 1324 | ! iwmax=nlev_ec+1 |
---|
| 1325 | ! do 60 i=0,nxmin1 |
---|
| 1326 | ! do 60 j=0,nymin1 |
---|
| 1327 | !60 wwh(i,j,nlev_ec+1)=0. |
---|
| 1328 | ! endif |
---|
| 1329 | |
---|
| 1330 | ! For global fields, assign the leftmost data column also to the rightmost |
---|
| 1331 | ! data column; if required, shift whole grid by nxshift grid points |
---|
| 1332 | ! |
---|
| 1333 | ! FLEXPART_WRF - all "global" stuff is turned off |
---|
| 1334 | !************************************************************************* |
---|
| 1335 | |
---|
| 1336 | ! if (xglobal) then |
---|
| 1337 | ! call shift_field_0(ewss,nxfield,ny) |
---|
| 1338 | ! call shift_field_0(nsss,nxfield,ny) |
---|
| 1339 | ! call shift_field_0(oro,nxfield,ny) |
---|
| 1340 | ! call shift_field_0(excessoro,nxfield,ny) |
---|
| 1341 | ! call shift_field_0(lsm,nxfield,ny) |
---|
| 1342 | ! call shift_field(ps,nxfield,ny,1,1,2,n) |
---|
| 1343 | ! call shift_field(sd,nxfield,ny,1,1,2,n) |
---|
| 1344 | ! call shift_field(msl,nxfield,ny,1,1,2,n) |
---|
| 1345 | ! call shift_field(tcc,nxfield,ny,1,1,2,n) |
---|
| 1346 | ! call shift_field(u10,nxfield,ny,1,1,2,n) |
---|
| 1347 | ! call shift_field(v10,nxfield,ny,1,1,2,n) |
---|
| 1348 | ! call shift_field(tt2,nxfield,ny,1,1,2,n) |
---|
| 1349 | ! call shift_field(td2,nxfield,ny,1,1,2,n) |
---|
| 1350 | ! call shift_field(lsprec,nxfield,ny,1,1,2,n) |
---|
| 1351 | ! call shift_field(convprec,nxfield,ny,1,1,2,n) |
---|
| 1352 | ! call shift_field(sshf,nxfield,ny,1,1,2,n) |
---|
| 1353 | ! call shift_field(ssr,nxfield,ny,1,1,2,n) |
---|
| 1354 | ! call shift_field(tth,nxfield,ny,nuvzmax,nuvz,2,n) |
---|
| 1355 | ! call shift_field(qvh,nxfield,ny,nuvzmax,nuvz,2,n) |
---|
| 1356 | ! call shift_field(uuh,nxfield,ny,nuvzmax,nuvz,1,1) |
---|
| 1357 | ! call shift_field(vvh,nxfield,ny,nuvzmax,nuvz,1,1) |
---|
| 1358 | ! call shift_field(wwh,nxfield,ny,nwzmax,nwz,1,1) |
---|
| 1359 | ! endif |
---|
| 1360 | |
---|
| 1361 | ! CALCULATE SURFSTR |
---|
| 1362 | if(sfc_option .eq. sfc_option_diagnosed) then |
---|
| 1363 | do i=0,nxmin1 |
---|
| 1364 | do j=0,nymin1 |
---|
| 1365 | surfstr(i,j,1,n)=sqrt(ewss(i,j)**2+nsss(i,j)**2) |
---|
| 1366 | enddo |
---|
| 1367 | enddo |
---|
| 1368 | strswitch=.false. |
---|
| 1369 | endif |
---|
| 1370 | |
---|
| 1371 | if ((.not.hflswitch).or.(.not.strswitch)) then |
---|
| 1372 | write(*,*) 'WARNING: No (or incomplete) flux data ' // & |
---|
| 1373 | 'contained in WRF output file ', & |
---|
| 1374 | wfname(indj) |
---|
| 1375 | |
---|
| 1376 | ! CALCULATE USTAR AND SSHF USING THE PROFILE METHOD |
---|
| 1377 | ! As ECMWF has increased the model resolution, such that now the first model |
---|
| 1378 | ! level is at about 10 m (where 10-m wind is given), use the 2nd ECMWF level |
---|
| 1379 | ! (3rd model level in FLEXPART) for the profile method |
---|
| 1380 | ! |
---|
| 1381 | ! FLEXPART_WRF - use k=(2+add_sfc_level) here instead of k=3 |
---|
| 1382 | !*************************************************************************** |
---|
| 1383 | |
---|
| 1384 | k = 2 + add_sfc_level |
---|
| 1385 | do i=0,nxmin1 |
---|
| 1386 | do j=0,nymin1 |
---|
| 1387 | ! plev1=akz(3)+bkz(3)*ps(i,j,1,n) |
---|
| 1388 | plev1=pph(i,j,k,n) |
---|
| 1389 | pmean=0.5*(ps(i,j,1,n)+plev1) |
---|
| 1390 | tv=tth(i,j,k,n)*(1.+0.61*qvh(i,j,k,n)) |
---|
| 1391 | fu=-r_air*tv/ga/pmean |
---|
| 1392 | hlev1=fu*(plev1-ps(i,j,1,n)) ! HEIGTH OF FIRST MODEL LAYER |
---|
| 1393 | ff10m= sqrt(u10(i,j,1,n)**2+v10(i,j,1,n)**2) |
---|
| 1394 | fflev1=sqrt(uuh(i,j,k)**2+vvh(i,j,k)**2) |
---|
| 1395 | call pbl_profile(ps(i,j,1,n),td2(i,j,1,n),hlev1, & |
---|
| 1396 | tt2(i,j,1,n),tth(i,j,k,n),ff10m,fflev1, & |
---|
| 1397 | surfstr(i,j,1,n),sshf(i,j,1,n)) |
---|
| 1398 | if(sshf(i,j,1,n).gt.200.) sshf(i,j,1,n)=200. |
---|
| 1399 | if(sshf(i,j,1,n).lt.-400.) sshf(i,j,1,n)=-400. |
---|
| 1400 | enddo |
---|
| 1401 | enddo |
---|
| 1402 | endif |
---|
| 1403 | |
---|
| 1404 | |
---|
| 1405 | ! Assign 10 m wind to model level at eta=1.0 to have one additional model |
---|
| 1406 | ! level at the ground |
---|
| 1407 | ! Specific humidity is taken the same as at one level above |
---|
| 1408 | ! Temperature is taken as 2 m temperature |
---|
| 1409 | ! |
---|
| 1410 | ! Note that the uuh, vvh, tth, & qvh data have already been shifted |
---|
| 1411 | ! upwards by one level, when they were read in. |
---|
| 1412 | !************************************************************************** |
---|
| 1413 | |
---|
| 1414 | if (add_sfc_level .eq. 1) then |
---|
| 1415 | do j = 0, nymin1 |
---|
| 1416 | do i = 0, nxmin1 |
---|
| 1417 | uuh(i,j,1) = u10(i,j,1,n) |
---|
| 1418 | vvh(i,j,1) = v10(i,j,1,n) |
---|
| 1419 | tth(i,j,1,n) = tt2(i,j,1,n) |
---|
| 1420 | ptth(i,j,1,n) = ptth(i,j,2,n) |
---|
| 1421 | qvh(i,j,1,n) = qvh(i,j,2,n) |
---|
| 1422 | tkeh(i,j,1,n)=tkeh(i,j,2,n) |
---|
| 1423 | ! pressure at 2 m AGL |
---|
| 1424 | pph(i,j,1,n) = 0.99976*ps(i,j,1,n) |
---|
| 1425 | ! height (MSL) at ground level (shift it down) |
---|
| 1426 | zzh(i,j,1,n) = zzh(i,j,2,n) |
---|
| 1427 | ! height (MSL) at top of the added level |
---|
| 1428 | zzh(i,j,2,n) = zzh(i,j,1,n) + 4.0 |
---|
| 1429 | if (hmix(i,j,1,n).lt.hmixmin) hmix(i,j,1,n)=hmixmin |
---|
| 1430 | |
---|
| 1431 | enddo |
---|
| 1432 | enddo |
---|
| 1433 | end if |
---|
| 1434 | |
---|
| 1435 | |
---|
| 1436 | do i=0,nxmax-1 |
---|
| 1437 | do j=0,nymax-1 |
---|
| 1438 | do k=1,nuvzmax |
---|
| 1439 | u_wrf(i,j,k,n)=uuh(i,j,k) |
---|
| 1440 | v_wrf(i,j,k,n)=vvh(i,j,k) |
---|
| 1441 | enddo |
---|
| 1442 | enddo |
---|
| 1443 | enddo |
---|
| 1444 | |
---|
| 1445 | do i=0,nxmax-1 |
---|
| 1446 | do j=0,nymax-1 |
---|
| 1447 | do k=1,nwzmax |
---|
| 1448 | w_wrf(i,j,k,n)=wwh(i,j,k) |
---|
| 1449 | enddo |
---|
| 1450 | enddo |
---|
| 1451 | enddo |
---|
| 1452 | |
---|
| 1453 | |
---|
| 1454 | |
---|
| 1455 | |
---|
| 1456 | |
---|
| 1457 | return |
---|
| 1458 | end subroutine readwind |
---|
| 1459 | |
---|