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