[2128] | 1 | MODULE obs_grid |
---|
| 2 | !!====================================================================== |
---|
| 3 | !! *** MODULE obs_grid *** |
---|
| 4 | !! Observation diagnostics: Various tools for grid searching etc. |
---|
| 5 | !!====================================================================== |
---|
| 6 | !!---------------------------------------------------------------------- |
---|
| 7 | !! obs_grid_search : Find i,j on the ORCA grid from lat,lon |
---|
| 8 | !! obs_level_search : Find level from depth |
---|
| 9 | !! obs_zlevel_search : Find depth level from observed depth |
---|
| 10 | !! obs_tlevel_search : Find temperature level from observed temp |
---|
| 11 | !! obs_rlevel_search : Find density level from observed density |
---|
| 12 | !!---------------------------------------------------------------------- |
---|
| 13 | !! * Modules used |
---|
| 14 | USE par_kind, ONLY : & ! Precision variables |
---|
| 15 | & wp |
---|
| 16 | USE par_oce, ONLY : & ! Ocean parameters |
---|
| 17 | & jpk, & |
---|
| 18 | & jpni, & |
---|
| 19 | & jpnj, & |
---|
| 20 | & jpnij |
---|
| 21 | USE dom_oce ! Ocean space and time domain variables |
---|
| 22 | USE obs_mpp, ONLY : & ! MPP support routines for observation diagnostics |
---|
| 23 | & obs_mpp_find_obs_proc, & |
---|
| 24 | & mpp_global_max, & |
---|
| 25 | & obs_mpp_max_integer |
---|
| 26 | USE phycst, ONLY : & ! Physical constants |
---|
| 27 | & rad |
---|
| 28 | USE obs_utils, ONLY : & ! Observation operator utility functions |
---|
| 29 | & grt_cir_dis, & |
---|
| 30 | & chkerr |
---|
| 31 | USE in_out_manager ! Printing support |
---|
| 32 | USE netcdf |
---|
| 33 | USE obs_const, ONLY : & |
---|
| 34 | & obfillflt ! Fillvalue |
---|
[2715] | 35 | USE lib_mpp, ONLY : & |
---|
| 36 | & ctl_warn, ctl_stop |
---|
| 37 | |
---|
[2128] | 38 | IMPLICIT NONE |
---|
| 39 | |
---|
| 40 | !! * Routine accessibility |
---|
| 41 | PUBLIC obs_grid_setup, & ! Setup grid searching |
---|
| 42 | & obs_grid_search, & ! Find i, j on the ORCA grid from lat, lon |
---|
| 43 | & obs_grid_deallocate, & ! Deallocate the look up table |
---|
| 44 | & obs_level_search ! Find level from depth |
---|
| 45 | |
---|
| 46 | PRIVATE linquad, & ! Determine whether a point lies within a cell |
---|
| 47 | & maxdist, & ! Find the maximum distance between 2 pts in a cell |
---|
[2358] | 48 | & obs_grd_bruteforce, & ! Find i, j on the ORCA grid from lat, lon |
---|
| 49 | & obs_grd_lookup ! Find i, j on the ORCA grid from lat, lon quicker |
---|
[2128] | 50 | |
---|
| 51 | !!* Module variables |
---|
| 52 | |
---|
| 53 | !! Default values |
---|
[6140] | 54 | REAL, PUBLIC :: rn_gridsearchres = 0.5 ! Resolution of grid |
---|
[2128] | 55 | INTEGER, PRIVATE :: gsearch_nlons_def ! Num of longitudes |
---|
| 56 | INTEGER, PRIVATE :: gsearch_nlats_def ! Num of latitudes |
---|
| 57 | REAL(wp), PRIVATE :: gsearch_lonmin_def ! Min longitude |
---|
| 58 | REAL(wp), PRIVATE :: gsearch_latmin_def ! Min latitude |
---|
| 59 | REAL(wp), PRIVATE :: gsearch_dlon_def ! Lon spacing |
---|
| 60 | REAL(wp), PRIVATE :: gsearch_dlat_def ! Lat spacing |
---|
| 61 | !! Variable versions |
---|
| 62 | INTEGER, PRIVATE :: nlons ! Num of longitudes |
---|
| 63 | INTEGER, PRIVATE :: nlats ! Num of latitudes |
---|
| 64 | REAL(wp), PRIVATE :: lonmin ! Min longitude |
---|
| 65 | REAL(wp), PRIVATE :: latmin ! Min latitude |
---|
| 66 | REAL(wp), PRIVATE :: dlon ! Lon spacing |
---|
| 67 | REAL(wp), PRIVATE :: dlat ! Lat spacing |
---|
| 68 | |
---|
| 69 | INTEGER, PRIVATE :: maxxdiff, maxydiff ! Max diffs between model points |
---|
| 70 | INTEGER, PRIVATE :: limxdiff, limydiff |
---|
| 71 | |
---|
| 72 | ! Data storage |
---|
| 73 | REAL(wp), PRIVATE, DIMENSION(:,:), ALLOCATABLE :: & |
---|
| 74 | & lons, & |
---|
| 75 | & lats |
---|
| 76 | INTEGER, PRIVATE, DIMENSION(:,:), ALLOCATABLE :: & |
---|
| 77 | & ixpos, & |
---|
| 78 | & iypos, & |
---|
[2358] | 79 | & iprocn |
---|
[2128] | 80 | |
---|
| 81 | ! Switches |
---|
| 82 | LOGICAL, PUBLIC :: ln_grid_search_lookup ! Use lookup table to speed up grid search |
---|
| 83 | LOGICAL, PUBLIC :: ln_grid_global ! Use global distribution of observations |
---|
| 84 | CHARACTER(LEN=44), PUBLIC :: & |
---|
[6140] | 85 | & cn_gridsearchfile ! file name head for grid search lookup |
---|
[2128] | 86 | |
---|
[2287] | 87 | !!---------------------------------------------------------------------- |
---|
[9598] | 88 | !! NEMO/OCE 4.0 , NEMO Consortium (2018) |
---|
[2287] | 89 | !! $Id$ |
---|
[10068] | 90 | !! Software governed by the CeCILL license (see ./LICENSE) |
---|
[2287] | 91 | !!---------------------------------------------------------------------- |
---|
| 92 | |
---|
[2128] | 93 | CONTAINS |
---|
| 94 | |
---|
| 95 | SUBROUTINE obs_grid_search( kobsin, plam, pphi, kobsi, kobsj, kproc, & |
---|
| 96 | & cdgrid ) |
---|
| 97 | !!---------------------------------------------------------------------- |
---|
| 98 | !! *** ROUTINE obs_grid_search *** |
---|
| 99 | !! |
---|
| 100 | !! ** Purpose : Search local gridpoints to find the grid box containing |
---|
| 101 | !! the observations calls either |
---|
[2358] | 102 | !! obs_grd_bruteforce - the original brute force search |
---|
[2128] | 103 | !! or |
---|
[2358] | 104 | !! obs_grd_lookup - uses a lookup table to do a fast |
---|
[2128] | 105 | !!search |
---|
| 106 | !!History : |
---|
| 107 | !! ! 2007-12 (D. Lea) |
---|
| 108 | !!------------------------------------------------------------------------ |
---|
| 109 | |
---|
| 110 | !! * Arguments |
---|
| 111 | INTEGER :: & |
---|
| 112 | & kobsin ! Size of the observation arrays |
---|
| 113 | REAL(KIND=wp), DIMENSION(kobsin), INTENT(IN) :: & |
---|
| 114 | & plam, & ! Longitude of obsrvations |
---|
| 115 | & pphi ! Latitude of observations |
---|
| 116 | INTEGER, DIMENSION(kobsin), INTENT(OUT) :: & |
---|
| 117 | & kobsi, & ! I-index of observations |
---|
| 118 | & kobsj, & ! J-index of observations |
---|
| 119 | & kproc ! Processor number of observations |
---|
| 120 | CHARACTER(LEN=1) :: & |
---|
| 121 | & cdgrid ! Grid to search |
---|
| 122 | |
---|
| 123 | IF(kobsin > 0) THEN |
---|
| 124 | |
---|
| 125 | IF ( ln_grid_search_lookup .AND. ( cdgrid == 'T' ) ) THEN |
---|
[2358] | 126 | CALL obs_grd_lookup( kobsin, plam, pphi, & |
---|
[2128] | 127 | & kobsi, kobsj, kproc ) |
---|
| 128 | ELSE |
---|
| 129 | IF ( cdgrid == 'T' ) THEN |
---|
[2358] | 130 | CALL obs_grd_bruteforce( jpi, jpj, jpiglo, jpjglo, & |
---|
[9523] | 131 | & 1, nlci, 1, nlcj, & |
---|
[2128] | 132 | & nproc, jpnij, & |
---|
| 133 | & glamt, gphit, tmask, & |
---|
| 134 | & kobsin, plam, pphi, & |
---|
| 135 | & kobsi, kobsj, kproc ) |
---|
| 136 | ELSEIF ( cdgrid == 'U' ) THEN |
---|
[2358] | 137 | CALL obs_grd_bruteforce( jpi, jpj, jpiglo, jpjglo, & |
---|
[9523] | 138 | & 1, nlci, 1, nlcj, & |
---|
[2128] | 139 | & nproc, jpnij, & |
---|
| 140 | & glamu, gphiu, umask, & |
---|
| 141 | & kobsin, plam, pphi, & |
---|
| 142 | & kobsi, kobsj, kproc ) |
---|
| 143 | ELSEIF ( cdgrid == 'V' ) THEN |
---|
[2358] | 144 | CALL obs_grd_bruteforce( jpi, jpj, jpiglo, jpjglo, & |
---|
[9523] | 145 | & 1, nlci, 1, nlcj, & |
---|
[2128] | 146 | & nproc, jpnij, & |
---|
| 147 | & glamv, gphiv, vmask, & |
---|
| 148 | & kobsin, plam, pphi, & |
---|
| 149 | & kobsi, kobsj, kproc ) |
---|
| 150 | ELSEIF ( cdgrid == 'F' ) THEN |
---|
[2358] | 151 | CALL obs_grd_bruteforce( jpi, jpj, jpiglo, jpjglo, & |
---|
[9523] | 152 | & 1, nlci, 1, nlcj, & |
---|
[2128] | 153 | & nproc, jpnij, & |
---|
| 154 | & glamf, gphif, fmask, & |
---|
| 155 | & kobsin, plam, pphi, & |
---|
| 156 | & kobsi, kobsj, kproc ) |
---|
| 157 | ELSE |
---|
| 158 | CALL ctl_stop( 'Grid not supported' ) |
---|
| 159 | ENDIF |
---|
| 160 | ENDIF |
---|
| 161 | |
---|
| 162 | ENDIF |
---|
| 163 | |
---|
| 164 | END SUBROUTINE obs_grid_search |
---|
| 165 | |
---|
[2358] | 166 | #include "obs_grd_bruteforce.h90" |
---|
[2128] | 167 | |
---|
[2358] | 168 | SUBROUTINE obs_grd_lookup( kobs, plam, pphi, kobsi, kobsj, kproc ) |
---|
[2128] | 169 | !!---------------------------------------------------------------------- |
---|
[2358] | 170 | !! *** ROUTINE obs_grid_lookup *** |
---|
[2128] | 171 | !! |
---|
| 172 | !! ** Purpose : Search local gridpoints to find the grid box containing |
---|
[2358] | 173 | !! the observations (much faster then obs_grd_bruteforce) |
---|
[2128] | 174 | !! |
---|
| 175 | !! ** Method : Call to linquad |
---|
| 176 | !! |
---|
| 177 | !! ** Action : Return kproc holding the observation and kiobsi,kobsj |
---|
| 178 | !! valid on kproc=nproc processor only. |
---|
| 179 | !! |
---|
| 180 | !! History : |
---|
| 181 | !! ! 2007-12 (D. Lea) new routine based on obs_grid_search |
---|
| 182 | !!! updated with fixes from new version of obs_grid_search_bruteforce |
---|
| 183 | !!! speeded up where points are not near a "difficult" region like an edge |
---|
| 184 | !!---------------------------------------------------------------------- |
---|
| 185 | |
---|
| 186 | !! * Arguments |
---|
| 187 | INTEGER :: kobs ! Size of the observation arrays |
---|
| 188 | REAL(KIND=wp), DIMENSION(kobs), INTENT(IN) :: & |
---|
| 189 | & plam, & ! Longitude of obsrvations |
---|
| 190 | & pphi ! Latitude of observations |
---|
| 191 | INTEGER, DIMENSION(kobs), INTENT(OUT) :: & |
---|
| 192 | & kobsi, & ! I-index of observations |
---|
| 193 | & kobsj, & ! J-index of observations |
---|
| 194 | & kproc ! Processor number of observations |
---|
| 195 | |
---|
| 196 | !! * Local declarations |
---|
| 197 | REAL(KIND=wp), DIMENSION(:), ALLOCATABLE :: & |
---|
| 198 | & zplam |
---|
| 199 | REAL(wp) :: zlammax |
---|
| 200 | REAL(wp) :: zlam |
---|
| 201 | INTEGER :: ji |
---|
| 202 | INTEGER :: jj |
---|
| 203 | INTEGER :: jk |
---|
| 204 | INTEGER :: jo |
---|
| 205 | INTEGER :: isx |
---|
| 206 | INTEGER :: isy |
---|
| 207 | INTEGER :: jimin |
---|
| 208 | INTEGER :: jimax |
---|
| 209 | INTEGER :: jjmin |
---|
| 210 | INTEGER :: jjmax |
---|
| 211 | INTEGER :: jojimin |
---|
| 212 | INTEGER :: jojimax |
---|
| 213 | INTEGER :: jojjmin |
---|
| 214 | INTEGER :: jojjmax |
---|
| 215 | INTEGER :: ipx1 |
---|
| 216 | INTEGER :: ipy1 |
---|
| 217 | INTEGER :: ip |
---|
| 218 | INTEGER :: jp |
---|
| 219 | INTEGER :: ipx |
---|
| 220 | INTEGER :: ipy |
---|
| 221 | INTEGER :: ipmx |
---|
| 222 | INTEGER :: jlon |
---|
| 223 | INTEGER :: jlat |
---|
| 224 | INTEGER :: joffset |
---|
| 225 | INTEGER :: jostride |
---|
| 226 | REAL(KIND=wp), DIMENSION(:,:), ALLOCATABLE :: & |
---|
| 227 | & zlamg, & |
---|
| 228 | & zphig, & |
---|
| 229 | & zmskg, & |
---|
| 230 | & zphitmax,& |
---|
| 231 | & zphitmin,& |
---|
| 232 | & zlamtmax,& |
---|
| 233 | & zlamtmin |
---|
| 234 | LOGICAL, DIMENSION(:,:), ALLOCATABLE :: & |
---|
| 235 | & llinvalidcell |
---|
| 236 | REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & |
---|
| 237 | & zlamtm, & |
---|
| 238 | & zphitm |
---|
| 239 | LOGICAL :: llfourflag |
---|
| 240 | INTEGER :: ifourflagcountt |
---|
| 241 | INTEGER :: ifourflagcountf |
---|
| 242 | INTEGER, DIMENSION(5) :: ifourflagcountr |
---|
| 243 | |
---|
| 244 | !----------------------------------------------------------------------- |
---|
| 245 | ! Define grid for grid search |
---|
| 246 | !----------------------------------------------------------------------- |
---|
| 247 | IF (ln_grid_global) THEN |
---|
| 248 | jlon = jpiglo |
---|
| 249 | jlat = jpjglo |
---|
| 250 | joffset = nproc |
---|
| 251 | jostride = jpnij |
---|
| 252 | ELSE |
---|
| 253 | jlon = jpi |
---|
| 254 | jlat = jpj |
---|
| 255 | joffset = 0 |
---|
| 256 | jostride = 1 |
---|
| 257 | ENDIF |
---|
| 258 | !----------------------------------------------------------------------- |
---|
| 259 | ! Set up data for grid search |
---|
| 260 | !----------------------------------------------------------------------- |
---|
| 261 | ALLOCATE( & |
---|
| 262 | & zlamg(jlon,jlat), & |
---|
| 263 | & zphig(jlon,jlat), & |
---|
| 264 | & zmskg(jlon,jlat), & |
---|
| 265 | & zphitmax(jlon-1,jlat-1), & |
---|
| 266 | & zphitmin(jlon-1,jlat-1), & |
---|
| 267 | & zlamtmax(jlon-1,jlat-1), & |
---|
| 268 | & zlamtmin(jlon-1,jlat-1), & |
---|
| 269 | & llinvalidcell(jlon-1,jlat-1), & |
---|
| 270 | & zlamtm(4,jlon-1,jlat-1), & |
---|
| 271 | & zphitm(4,jlon-1,jlat-1) & |
---|
| 272 | & ) |
---|
| 273 | !----------------------------------------------------------------------- |
---|
| 274 | ! Copy data to local arrays |
---|
| 275 | !----------------------------------------------------------------------- |
---|
| 276 | IF (ln_grid_global) THEN |
---|
| 277 | zlamg(:,:) = -1.e+10 |
---|
| 278 | zphig(:,:) = -1.e+10 |
---|
| 279 | zmskg(:,:) = -1.e+10 |
---|
| 280 | ! Add various grids here. |
---|
[9523] | 281 | DO jj = 1, nlcj |
---|
| 282 | DO ji = 1, nlci |
---|
[2128] | 283 | zlamg(mig(ji),mjg(jj)) = glamt(ji,jj) |
---|
| 284 | zphig(mig(ji),mjg(jj)) = gphit(ji,jj) |
---|
| 285 | zmskg(mig(ji),mjg(jj)) = tmask(ji,jj,1) |
---|
| 286 | END DO |
---|
| 287 | END DO |
---|
| 288 | CALL mpp_global_max( zlamg ) |
---|
| 289 | CALL mpp_global_max( zphig ) |
---|
| 290 | CALL mpp_global_max( zmskg ) |
---|
| 291 | ELSE |
---|
| 292 | ! Add various grids here. |
---|
| 293 | DO jj = 1, jlat |
---|
| 294 | DO ji = 1, jlon |
---|
| 295 | zlamg(ji,jj) = glamt(ji,jj) |
---|
| 296 | zphig(ji,jj) = gphit(ji,jj) |
---|
| 297 | zmskg(ji,jj) = tmask(ji,jj,1) |
---|
| 298 | END DO |
---|
| 299 | END DO |
---|
| 300 | ENDIF |
---|
| 301 | !----------------------------------------------------------------------- |
---|
| 302 | ! Copy longitudes |
---|
| 303 | !----------------------------------------------------------------------- |
---|
| 304 | ALLOCATE( & |
---|
| 305 | & zplam(kobs) & |
---|
| 306 | & ) |
---|
| 307 | DO jo = 1, kobs |
---|
| 308 | zplam(jo) = plam(jo) |
---|
| 309 | END DO |
---|
| 310 | !----------------------------------------------------------------------- |
---|
| 311 | ! Set default values for output |
---|
| 312 | !----------------------------------------------------------------------- |
---|
| 313 | kproc(:) = -1 |
---|
| 314 | kobsi(:) = -1 |
---|
| 315 | kobsj(:) = -1 |
---|
| 316 | !----------------------------------------------------------------------- |
---|
| 317 | ! Copy grid positions to temporary arrays and renormalize to 0 to 360. |
---|
| 318 | !----------------------------------------------------------------------- |
---|
| 319 | DO jj = 1, jlat-1 |
---|
| 320 | DO ji = 1, jlon-1 |
---|
| 321 | zlamtm(1,ji,jj) = zlamg(ji ,jj ) |
---|
| 322 | zphitm(1,ji,jj) = zphig(ji ,jj ) |
---|
| 323 | zlamtm(2,ji,jj) = zlamg(ji+1,jj ) |
---|
| 324 | zphitm(2,ji,jj) = zphig(ji+1,jj ) |
---|
| 325 | zlamtm(3,ji,jj) = zlamg(ji+1,jj+1) |
---|
| 326 | zphitm(3,ji,jj) = zphig(ji+1,jj+1) |
---|
| 327 | zlamtm(4,ji,jj) = zlamg(ji ,jj+1) |
---|
| 328 | zphitm(4,ji,jj) = zphig(ji ,jj+1) |
---|
| 329 | END DO |
---|
| 330 | END DO |
---|
| 331 | WHERE ( zlamtm(:,:,:) < 0.0_wp ) |
---|
| 332 | zlamtm(:,:,:) = zlamtm(:,:,:) + 360.0_wp |
---|
| 333 | END WHERE |
---|
| 334 | WHERE ( zlamtm(:,:,:) > 360.0_wp ) |
---|
| 335 | zlamtm(:,:,:) = zlamtm(:,:,:) - 360.0_wp |
---|
| 336 | END WHERE |
---|
| 337 | !----------------------------------------------------------------------- |
---|
| 338 | ! Handle case of the wraparound; beware, not working with orca180 |
---|
| 339 | !----------------------------------------------------------------------- |
---|
| 340 | DO jj = 1, jlat-1 |
---|
| 341 | DO ji = 1, jlon-1 |
---|
| 342 | zlammax = MAXVAL( zlamtm(:,ji,jj) ) |
---|
| 343 | WHERE (zlammax - zlamtm(:, ji, jj) > 180 ) & |
---|
| 344 | & zlamtm(:,ji,jj) = zlamtm(:,ji,jj) + 360._wp |
---|
| 345 | zphitmax(ji,jj) = MAXVAL(zphitm(:,ji,jj)) |
---|
| 346 | zphitmin(ji,jj) = MINVAL(zphitm(:,ji,jj)) |
---|
| 347 | zlamtmax(ji,jj) = MAXVAL(zlamtm(:,ji,jj)) |
---|
| 348 | zlamtmin(ji,jj) = MINVAL(zlamtm(:,ji,jj)) |
---|
| 349 | END DO |
---|
| 350 | END DO |
---|
| 351 | !----------------------------------------------------------------------- |
---|
| 352 | ! Search for boxes with only land points mark them invalid |
---|
| 353 | !----------------------------------------------------------------------- |
---|
| 354 | llinvalidcell(:,:) = .FALSE. |
---|
| 355 | DO jj = 1, jlat-1 |
---|
| 356 | DO ji = 1, jlon-1 |
---|
| 357 | llinvalidcell(ji,jj) = & |
---|
| 358 | & zmskg(ji ,jj ) == 0.0_wp .AND. & |
---|
| 359 | & zmskg(ji+1,jj ) == 0.0_wp .AND. & |
---|
| 360 | & zmskg(ji+1,jj+1) == 0.0_wp .AND. & |
---|
| 361 | & zmskg(ji ,jj+1) == 0.0_wp |
---|
| 362 | END DO |
---|
| 363 | END DO |
---|
| 364 | |
---|
[2358] | 365 | if(lwp) WRITE(numout,*) 'obs_grid_lookup do coordinate search using lookup table' |
---|
[2128] | 366 | |
---|
| 367 | !----------------------------------------------------------------------- |
---|
| 368 | ! Do coordinate search using lookup table with local searches. |
---|
| 369 | ! - For land points kproc is set to number of the processor + 1000000 |
---|
| 370 | ! and we continue the search. |
---|
| 371 | ! - For ocean points kproc is set to the number of the processor |
---|
| 372 | ! and we stop the search. |
---|
| 373 | !----------------------------------------------------------------------- |
---|
| 374 | ifourflagcountt = 0 |
---|
| 375 | ifourflagcountf = 0 |
---|
| 376 | ifourflagcountr(:) = 0 |
---|
| 377 | |
---|
| 378 | !------------------------------------------------------------------------ |
---|
| 379 | ! Master loop for grid search |
---|
| 380 | !------------------------------------------------------------------------ |
---|
| 381 | |
---|
| 382 | gpkobs: DO jo = 1+joffset, kobs, jostride |
---|
| 383 | ! Normal case |
---|
| 384 | ! specify 4 points which surround the lat lon of interest |
---|
| 385 | ! x i,j+1 x i+1, j+1 |
---|
| 386 | ! |
---|
| 387 | ! |
---|
| 388 | ! * lon,lat |
---|
| 389 | ! x i,j x i+1,j |
---|
| 390 | |
---|
| 391 | ! bottom corner point |
---|
| 392 | ipx1 = INT( ( zplam(jo) - lonmin ) / dlon + 1.0 ) |
---|
| 393 | ipy1 = INT( ( pphi (jo) - latmin ) / dlat + 1.0 ) |
---|
| 394 | |
---|
| 395 | ipx = ipx1 + 1 |
---|
| 396 | ipy = ipy1 + 1 |
---|
| 397 | |
---|
| 398 | ! flag for searching around four points separately |
---|
| 399 | ! default to false |
---|
| 400 | llfourflag = .FALSE. |
---|
| 401 | |
---|
| 402 | ! check for point fully outside of region |
---|
| 403 | IF ( (ipx1 > nlons) .OR. (ipy1 > nlats) .OR. & |
---|
| 404 | & (ipx < 1) .OR. (ipy < 1) ) THEN |
---|
| 405 | CYCLE |
---|
| 406 | ENDIF |
---|
| 407 | ! check wrap around |
---|
| 408 | IF ( (ipx > nlons) .OR. (ipy > nlats) .OR. & |
---|
| 409 | & (ipx1 < 1) .OR. (ipy1 < 1) ) THEN |
---|
| 410 | llfourflag=.TRUE. |
---|
| 411 | ifourflagcountr(1) = ifourflagcountr(1) + 1 |
---|
| 412 | ENDIF |
---|
| 413 | |
---|
[4990] | 414 | IF (.NOT. llfourflag) THEN |
---|
| 415 | IF (MAXVAL(ixpos(ipx1:ipx,ipy1:ipy)) == -1) CYCLE! cycle if no lookup points found |
---|
| 416 | ENDIF |
---|
[2128] | 417 | |
---|
| 418 | jimin = 0 |
---|
| 419 | jimax = 0 |
---|
| 420 | jjmin = 0 |
---|
| 421 | jjmax = 0 |
---|
| 422 | |
---|
| 423 | IF (.NOT. llfourflag) THEN |
---|
| 424 | |
---|
| 425 | ! calculate points range |
---|
| 426 | ! define a square region encompassing the four corner points |
---|
| 427 | ! do I need the -1 points? |
---|
| 428 | |
---|
| 429 | jojimin = MINVAL(ixpos(ipx1:ipx,ipy1:ipy)) - 1 |
---|
| 430 | jojimax = MAXVAL(ixpos(ipx1:ipx,ipy1:ipy)) + 1 |
---|
| 431 | jojjmin = MINVAL(iypos(ipx1:ipx,ipy1:ipy)) - 1 |
---|
| 432 | jojjmax = MAXVAL(iypos(ipx1:ipx,ipy1:ipy)) + 1 |
---|
| 433 | |
---|
| 434 | jimin = jojimin - 1 |
---|
| 435 | jimax = jojimax + 1 |
---|
| 436 | jjmin = jojjmin - 1 |
---|
| 437 | jjmax = jojjmax + 1 |
---|
| 438 | |
---|
| 439 | IF ( jojimin < 0 .OR. jojjmin < 0) THEN |
---|
| 440 | llfourflag = .TRUE. |
---|
| 441 | ifourflagcountr(2) = ifourflagcountr(2) + 1 |
---|
| 442 | ENDIF |
---|
| 443 | IF ( jojimax - jojimin > maxxdiff) THEN |
---|
| 444 | llfourflag = .TRUE. |
---|
| 445 | ifourflagcountr(3) = ifourflagcountr(3) + 1 |
---|
| 446 | ENDIF |
---|
| 447 | IF ( jojjmax - jojjmin > maxydiff) THEN |
---|
| 448 | llfourflag = .TRUE. |
---|
| 449 | ifourflagcountr(4) = ifourflagcountr(4) + 1 |
---|
| 450 | ENDIF |
---|
| 451 | |
---|
| 452 | ENDIF |
---|
| 453 | |
---|
| 454 | ipmx = 0 |
---|
| 455 | IF (llfourflag) ipmx = 1 |
---|
| 456 | |
---|
| 457 | IF (llfourflag) THEN |
---|
| 458 | ifourflagcountt = ifourflagcountt + 1 |
---|
| 459 | ELSE |
---|
| 460 | ifourflagcountf = ifourflagcountf + 1 |
---|
| 461 | ENDIF |
---|
| 462 | |
---|
| 463 | gridpointsn : DO ip = 0, ipmx |
---|
| 464 | DO jp = 0, ipmx |
---|
| 465 | |
---|
| 466 | IF ( kproc(jo) /= -1 ) EXIT gridpointsn |
---|
| 467 | |
---|
| 468 | ipx = ipx1 + ip |
---|
| 469 | ipy = ipy1 + jp |
---|
| 470 | |
---|
| 471 | IF (llfourflag) THEN |
---|
| 472 | |
---|
| 473 | ! deal with wrap around |
---|
| 474 | IF ( ipx > nlons ) ipx = 1 |
---|
| 475 | IF ( ipy > nlats ) ipy = 1 |
---|
| 476 | IF ( ipx < 1 ) ipx = nlons |
---|
| 477 | IF ( ipy < 1 ) ipy = nlats |
---|
| 478 | |
---|
| 479 | ! get i,j |
---|
| 480 | isx = ixpos(ipx,ipy) |
---|
| 481 | isy = iypos(ipx,ipy) |
---|
| 482 | |
---|
| 483 | ! estimate appropriate search region (use max/min values) |
---|
| 484 | jimin = isx - maxxdiff - 1 |
---|
| 485 | jimax = isx + maxxdiff + 1 |
---|
| 486 | jjmin = isy - maxydiff - 1 |
---|
| 487 | jjmax = isy + maxydiff + 1 |
---|
| 488 | |
---|
| 489 | ENDIF |
---|
| 490 | |
---|
| 491 | IF ( jimin < 1 ) jimin = 1 |
---|
| 492 | IF ( jimax > jlon-1 ) jimax = jlon-1 |
---|
| 493 | IF ( jjmin < 1 ) jjmin = 1 |
---|
| 494 | IF ( jjmax > jlat-1 ) jjmax = jlat-1 |
---|
| 495 | |
---|
| 496 | !--------------------------------------------------------------- |
---|
| 497 | ! Ensure that all observation longtiudes are between 0 and 360 |
---|
| 498 | !--------------------------------------------------------------- |
---|
| 499 | |
---|
| 500 | IF ( zplam(jo) < 0.0_wp ) zplam(jo) = zplam(jo) + 360.0_wp |
---|
| 501 | IF ( zplam(jo) > 360.0_wp ) zplam(jo) = zplam(jo) - 360.0_wp |
---|
| 502 | |
---|
| 503 | !--------------------------------------------------------------- |
---|
| 504 | ! Find observations which are on within 1e-6 of a grid point |
---|
| 505 | !--------------------------------------------------------------- |
---|
| 506 | |
---|
| 507 | gridloop: DO jj = jjmin, jjmax |
---|
| 508 | DO ji = jimin, jimax |
---|
| 509 | IF ( ABS( zphig(ji,jj) - pphi(jo) ) < 1e-6 ) THEN |
---|
| 510 | zlam = zlamg(ji,jj) |
---|
| 511 | IF ( zlam < 0.0_wp ) zlam = zlam + 360.0_wp |
---|
| 512 | IF ( zlam > 360.0_wp ) zlam = zlam - 360.0_wp |
---|
| 513 | IF ( ABS( zlam - zplam(jo) ) < 1e-6 ) THEN |
---|
| 514 | IF ( llinvalidcell(ji,jj) ) THEN |
---|
| 515 | kproc(jo) = nproc + 1000000 |
---|
| 516 | kobsi(jo) = ji + 1 |
---|
| 517 | kobsj(jo) = jj + 1 |
---|
| 518 | CYCLE |
---|
| 519 | ELSE |
---|
| 520 | kproc(jo) = nproc |
---|
| 521 | kobsi(jo) = ji + 1 |
---|
| 522 | kobsj(jo) = jj + 1 |
---|
| 523 | EXIT gridloop |
---|
| 524 | ENDIF |
---|
| 525 | ENDIF |
---|
| 526 | ENDIF |
---|
| 527 | END DO |
---|
| 528 | END DO gridloop |
---|
| 529 | |
---|
| 530 | !--------------------------------------------------------------- |
---|
| 531 | ! Ensure that all observation longtiudes are between -180/180 |
---|
| 532 | !--------------------------------------------------------------- |
---|
| 533 | |
---|
| 534 | IF ( zplam(jo) > 180 ) zplam(jo) = zplam(jo) - 360.0_wp |
---|
| 535 | |
---|
| 536 | IF ( kproc(jo) == -1 ) THEN |
---|
| 537 | |
---|
| 538 | ! Normal case |
---|
| 539 | gridpoints : DO jj = jjmin, jjmax |
---|
| 540 | DO ji = jimin, jimax |
---|
| 541 | |
---|
| 542 | |
---|
| 543 | IF ( ( zplam(jo) > zlamtmax(ji,jj) ) .OR. & |
---|
| 544 | & ( zplam(jo) < zlamtmin(ji,jj) ) ) CYCLE |
---|
| 545 | |
---|
| 546 | IF ( ABS( pphi(jo) ) < 85 ) THEN |
---|
| 547 | IF ( ( pphi(jo) > zphitmax(ji,jj) ) .OR. & |
---|
| 548 | & ( pphi(jo) < zphitmin(ji,jj) ) ) CYCLE |
---|
| 549 | ENDIF |
---|
| 550 | |
---|
| 551 | IF ( linquad( zplam(jo), pphi(jo), & |
---|
| 552 | & zlamtm(:,ji,jj), zphitm(:,ji,jj) ) ) THEN |
---|
| 553 | IF ( llinvalidcell(ji,jj) ) THEN |
---|
| 554 | kproc(jo) = nproc + 1000000 |
---|
| 555 | kobsi(jo) = ji + 1 |
---|
| 556 | kobsj(jo) = jj + 1 |
---|
| 557 | CYCLE |
---|
| 558 | ELSE |
---|
| 559 | kproc(jo) = nproc |
---|
| 560 | kobsi(jo) = ji + 1 |
---|
| 561 | kobsj(jo) = jj + 1 |
---|
| 562 | EXIT gridpoints |
---|
| 563 | ENDIF |
---|
| 564 | ENDIF |
---|
| 565 | |
---|
| 566 | END DO |
---|
| 567 | END DO gridpoints |
---|
| 568 | ENDIF |
---|
| 569 | |
---|
| 570 | ! In case of failure retry for obs. longtiude + 360. |
---|
| 571 | IF ( kproc(jo) == -1 ) THEN |
---|
| 572 | gridpoints_greenwich : DO jj = jjmin, jjmax |
---|
| 573 | DO ji = jimin, jimax |
---|
| 574 | |
---|
| 575 | IF ( ( zplam(jo)+360.0_wp > zlamtmax(ji,jj) ) .OR. & |
---|
| 576 | & ( zplam(jo)+360.0_wp < zlamtmin(ji,jj) ) ) CYCLE |
---|
| 577 | |
---|
| 578 | IF ( ABS( pphi(jo) ) < 85 ) THEN |
---|
| 579 | IF ( ( pphi(jo) > zphitmax(ji,jj) ) .OR. & |
---|
| 580 | & ( pphi(jo) < zphitmin(ji,jj) ) ) CYCLE |
---|
| 581 | ENDIF |
---|
| 582 | |
---|
| 583 | IF ( linquad( zplam(jo)+360.0_wp, pphi(jo), & |
---|
| 584 | & zlamtm(:,ji,jj), zphitm(:,ji,jj) ) ) THEN |
---|
| 585 | IF ( llinvalidcell(ji,jj) ) THEN |
---|
| 586 | kproc(jo) = nproc + 1000000 |
---|
| 587 | kobsi(jo) = ji + 1 |
---|
| 588 | kobsj(jo) = jj + 1 |
---|
| 589 | CYCLE |
---|
| 590 | ELSE |
---|
| 591 | kproc(jo) = nproc |
---|
| 592 | kobsi(jo) = ji + 1 |
---|
| 593 | kobsj(jo) = jj + 1 |
---|
| 594 | EXIT gridpoints_greenwich |
---|
| 595 | ENDIF |
---|
| 596 | ENDIF |
---|
| 597 | |
---|
| 598 | END DO |
---|
| 599 | END DO gridpoints_greenwich |
---|
| 600 | |
---|
| 601 | ENDIF ! kproc |
---|
| 602 | |
---|
| 603 | END DO |
---|
| 604 | END DO gridpointsn |
---|
| 605 | END DO gpkobs ! kobs |
---|
| 606 | |
---|
| 607 | !---------------------------------------------------------------------- |
---|
| 608 | ! Synchronize kproc on all processors |
---|
| 609 | !---------------------------------------------------------------------- |
---|
| 610 | IF ( ln_grid_global ) THEN |
---|
| 611 | CALL obs_mpp_max_integer( kproc, kobs ) |
---|
| 612 | CALL obs_mpp_max_integer( kobsi, kobs ) |
---|
| 613 | CALL obs_mpp_max_integer( kobsj, kobs ) |
---|
| 614 | ELSE |
---|
[6140] | 615 | CALL obs_mpp_find_obs_proc( kproc, kobs ) |
---|
[2128] | 616 | ENDIF |
---|
| 617 | |
---|
| 618 | WHERE( kproc(:) >= 1000000 ) |
---|
| 619 | kproc(:) = kproc(:) - 1000000 |
---|
| 620 | END WHERE |
---|
| 621 | |
---|
| 622 | DEALLOCATE( & |
---|
| 623 | & zlamg, & |
---|
| 624 | & zphig, & |
---|
| 625 | & zmskg, & |
---|
| 626 | & zphitmax, & |
---|
| 627 | & zphitmin, & |
---|
| 628 | & zlamtmax, & |
---|
| 629 | & zlamtmin, & |
---|
| 630 | & llinvalidcell, & |
---|
| 631 | & zlamtm, & |
---|
| 632 | & zphitm, & |
---|
| 633 | & zplam & |
---|
| 634 | & ) |
---|
| 635 | |
---|
[2358] | 636 | END SUBROUTINE obs_grd_lookup |
---|
[2128] | 637 | |
---|
| 638 | |
---|
| 639 | SUBROUTINE obs_grid_setup |
---|
| 640 | !!---------------------------------------------------------------------- |
---|
| 641 | !! *** ROUTINE obs_grid_setup *** |
---|
| 642 | !! |
---|
| 643 | !! ** Purpose : Setup a lookup table to reduce the searching required |
---|
| 644 | !! for converting lat lons to grid point location |
---|
| 645 | !! produces or reads in a preexisting file for use in |
---|
| 646 | !! obs_grid_search_lookup_local |
---|
| 647 | !! |
---|
| 648 | !! ** Method : calls obs_grid_search_bruteforce_local with a array |
---|
| 649 | !! of lats and lons |
---|
| 650 | !! |
---|
| 651 | !! History : |
---|
| 652 | !! ! 2007-12 (D. Lea) new routine |
---|
| 653 | !!---------------------------------------------------------------------- |
---|
| 654 | |
---|
| 655 | !! * Local declarations |
---|
| 656 | CHARACTER(LEN=15), PARAMETER :: & |
---|
| 657 | & cpname = 'obs_grid_setup' |
---|
| 658 | CHARACTER(LEN=40) :: cfname |
---|
| 659 | INTEGER :: ji |
---|
| 660 | INTEGER :: jj |
---|
| 661 | INTEGER :: jk |
---|
| 662 | INTEGER :: jo |
---|
| 663 | INTEGER :: idfile, idny, idnx, idxpos, idypos |
---|
| 664 | INTEGER :: idlat, idlon, fileexist |
---|
| 665 | INTEGER, DIMENSION(2) :: incdim |
---|
| 666 | CHARACTER(LEN=20) :: datestr=" ",timestr=" " |
---|
| 667 | REAL(wp) :: tmpx1, tmpx2, tmpy1, tmpy2 |
---|
| 668 | REAL(wp) :: meanxdiff, meanydiff |
---|
| 669 | REAL(wp) :: meanxdiff1, meanydiff1 |
---|
| 670 | REAL(wp) :: meanxdiff2, meanydiff2 |
---|
| 671 | INTEGER :: numx1, numx2, numy1, numy2, df |
---|
| 672 | INTEGER :: jimin, jimax, jjmin, jjmax |
---|
| 673 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: & |
---|
| 674 | & lonsi, & |
---|
| 675 | & latsi |
---|
| 676 | INTEGER, DIMENSION(:,:), ALLOCATABLE :: & |
---|
| 677 | & ixposi, & |
---|
| 678 | & iyposi, & |
---|
| 679 | & iproci |
---|
| 680 | INTEGER, PARAMETER :: histsize=90 |
---|
| 681 | INTEGER, DIMENSION(histsize) :: & |
---|
| 682 | & histx1, histx2, histy1, histy2 |
---|
| 683 | REAL, DIMENSION(histsize) :: & |
---|
| 684 | & fhistx1, fhistx2, fhisty1, fhisty2 |
---|
| 685 | REAL(wp) :: histtol |
---|
| 686 | |
---|
| 687 | IF (ln_grid_search_lookup) THEN |
---|
| 688 | |
---|
| 689 | WRITE(numout,*) 'Calling obs_grid_setup' |
---|
| 690 | |
---|
| 691 | IF(lwp) WRITE(numout,*) |
---|
[6140] | 692 | IF(lwp) WRITE(numout,*)'Grid search resolution : ', rn_gridsearchres |
---|
[2128] | 693 | |
---|
[6140] | 694 | gsearch_nlons_def = NINT( 360.0_wp / rn_gridsearchres ) |
---|
| 695 | gsearch_nlats_def = NINT( 180.0_wp / rn_gridsearchres ) |
---|
| 696 | gsearch_lonmin_def = -180.0_wp + 0.5_wp * rn_gridsearchres |
---|
| 697 | gsearch_latmin_def = -90.0_wp + 0.5_wp * rn_gridsearchres |
---|
| 698 | gsearch_dlon_def = rn_gridsearchres |
---|
| 699 | gsearch_dlat_def = rn_gridsearchres |
---|
[2128] | 700 | |
---|
| 701 | IF (lwp) THEN |
---|
| 702 | WRITE(numout,*)'Grid search gsearch_nlons_def = ',gsearch_nlons_def |
---|
| 703 | WRITE(numout,*)'Grid search gsearch_nlats_def = ',gsearch_nlats_def |
---|
| 704 | WRITE(numout,*)'Grid search gsearch_lonmin_def = ',gsearch_lonmin_def |
---|
| 705 | WRITE(numout,*)'Grid search gsearch_latmin_def = ',gsearch_latmin_def |
---|
| 706 | WRITE(numout,*)'Grid search gsearch_dlon_def = ',gsearch_dlon_def |
---|
| 707 | WRITE(numout,*)'Grid search gsearch_dlat_def = ',gsearch_dlat_def |
---|
| 708 | ENDIF |
---|
| 709 | |
---|
| 710 | IF ( ln_grid_global ) THEN |
---|
| 711 | WRITE(cfname, FMT="(A,'_',A)") & |
---|
[6140] | 712 | & TRIM(cn_gridsearchfile), 'global.nc' |
---|
[2128] | 713 | ELSE |
---|
| 714 | WRITE(cfname, FMT="(A,'_',I4.4,'of',I4.4,'by',I4.4,'.nc')") & |
---|
[6140] | 715 | & TRIM(cn_gridsearchfile), nproc, jpni, jpnj |
---|
[2128] | 716 | ENDIF |
---|
| 717 | |
---|
| 718 | fileexist=nf90_open( TRIM( cfname ), nf90_nowrite, & |
---|
| 719 | & idfile ) |
---|
| 720 | |
---|
| 721 | IF ( fileexist == nf90_noerr ) THEN |
---|
| 722 | |
---|
| 723 | ! read data |
---|
| 724 | ! initially assume size is as defined (to be fixed) |
---|
| 725 | |
---|
| 726 | WRITE(numout,*) 'Reading: ',cfname |
---|
| 727 | |
---|
| 728 | CALL chkerr( nf90_open( TRIM( cfname ), nf90_nowrite, idfile ), & |
---|
| 729 | & cpname, __LINE__ ) |
---|
| 730 | CALL chkerr( nf90_get_att( idfile, nf90_global, 'maxxdiff', maxxdiff ), & |
---|
| 731 | & cpname, __LINE__ ) |
---|
| 732 | CALL chkerr( nf90_get_att( idfile, nf90_global, 'maxydiff', maxydiff ), & |
---|
| 733 | & cpname, __LINE__ ) |
---|
| 734 | CALL chkerr( nf90_get_att( idfile, nf90_global, 'dlon', dlon ), & |
---|
| 735 | & cpname, __LINE__ ) |
---|
| 736 | CALL chkerr( nf90_get_att( idfile, nf90_global, 'dlat', dlat ), & |
---|
| 737 | & cpname, __LINE__ ) |
---|
| 738 | CALL chkerr( nf90_get_att( idfile, nf90_global, 'lonmin', lonmin ), & |
---|
| 739 | & cpname, __LINE__ ) |
---|
| 740 | CALL chkerr( nf90_get_att( idfile, nf90_global, 'latmin', latmin ), & |
---|
| 741 | & cpname, __LINE__ ) |
---|
| 742 | |
---|
| 743 | CALL chkerr( nf90_inq_dimid(idfile, 'nx' , idnx), & |
---|
| 744 | & cpname, __LINE__ ) |
---|
| 745 | CALL chkerr( nf90_inquire_dimension( idfile, idnx, len = nlons ), & |
---|
| 746 | & cpname, __LINE__ ) |
---|
| 747 | CALL chkerr( nf90_inq_dimid(idfile, 'ny' , idny), & |
---|
| 748 | & cpname, __LINE__ ) |
---|
| 749 | CALL chkerr( nf90_inquire_dimension( idfile, idny, len = nlats ), & |
---|
| 750 | & cpname, __LINE__ ) |
---|
| 751 | |
---|
| 752 | ALLOCATE( & |
---|
| 753 | & lons(nlons,nlats), & |
---|
| 754 | & lats(nlons,nlats), & |
---|
| 755 | & ixpos(nlons,nlats), & |
---|
| 756 | & iypos(nlons,nlats), & |
---|
[2358] | 757 | & iprocn(nlons,nlats) & |
---|
[2128] | 758 | & ) |
---|
| 759 | |
---|
| 760 | CALL chkerr( nf90_inq_varid( idfile, 'XPOS', idxpos ), & |
---|
| 761 | & cpname, __LINE__ ) |
---|
| 762 | CALL chkerr( nf90_get_var ( idfile, idxpos, ixpos), & |
---|
| 763 | & cpname, __LINE__ ) |
---|
| 764 | CALL chkerr( nf90_inq_varid( idfile, 'YPOS', idypos ), & |
---|
| 765 | & cpname, __LINE__ ) |
---|
| 766 | CALL chkerr( nf90_get_var ( idfile, idypos, iypos), & |
---|
| 767 | & cpname, __LINE__ ) |
---|
| 768 | |
---|
| 769 | CALL chkerr( nf90_close( idfile ), cpname, __LINE__ ) |
---|
| 770 | |
---|
| 771 | ! setup arrays |
---|
| 772 | |
---|
| 773 | DO ji = 1, nlons |
---|
| 774 | DO jj = 1, nlats |
---|
| 775 | lons(ji,jj) = lonmin + (ji-1) * dlon |
---|
| 776 | lats(ji,jj) = latmin + (jj-1) * dlat |
---|
| 777 | END DO |
---|
| 778 | END DO |
---|
| 779 | |
---|
| 780 | ! if we are not reading the file we need to create it |
---|
| 781 | ! create new obs grid search lookup file |
---|
| 782 | |
---|
| 783 | ELSE |
---|
| 784 | |
---|
| 785 | ! call obs_grid_search |
---|
| 786 | |
---|
| 787 | IF (lwp) THEN |
---|
| 788 | WRITE(numout,*) 'creating: ',cfname |
---|
| 789 | WRITE(numout,*) 'calling obs_grid_search: ',nlons*nlats |
---|
| 790 | ENDIF |
---|
| 791 | |
---|
| 792 | ! set parameters from default values |
---|
| 793 | nlons = gsearch_nlons_def |
---|
| 794 | nlats = gsearch_nlats_def |
---|
| 795 | lonmin = gsearch_lonmin_def |
---|
| 796 | latmin = gsearch_latmin_def |
---|
| 797 | dlon = gsearch_dlon_def |
---|
| 798 | dlat = gsearch_dlat_def |
---|
| 799 | |
---|
| 800 | ! setup arrays |
---|
| 801 | |
---|
| 802 | ALLOCATE( & |
---|
| 803 | & lonsi(nlons,nlats), & |
---|
| 804 | & latsi(nlons,nlats), & |
---|
| 805 | & ixposi(nlons,nlats), & |
---|
| 806 | & iyposi(nlons,nlats), & |
---|
| 807 | & iproci(nlons,nlats) & |
---|
| 808 | & ) |
---|
| 809 | |
---|
| 810 | DO ji = 1, nlons |
---|
| 811 | DO jj = 1, nlats |
---|
| 812 | lonsi(ji,jj) = lonmin + (ji-1) * dlon |
---|
| 813 | latsi(ji,jj) = latmin + (jj-1) * dlat |
---|
| 814 | END DO |
---|
| 815 | END DO |
---|
| 816 | |
---|
[2358] | 817 | CALL obs_grd_bruteforce( jpi, jpj, jpiglo, jpjglo, & |
---|
[9523] | 818 | & 1, nlci, 1, nlcj, & |
---|
[2358] | 819 | & nproc, jpnij, & |
---|
| 820 | & glamt, gphit, tmask, & |
---|
| 821 | & nlons*nlats, lonsi, latsi, & |
---|
| 822 | & ixposi, iyposi, iproci ) |
---|
[2128] | 823 | |
---|
| 824 | ! minimise file size by removing regions with no data from xypos file |
---|
| 825 | ! should be able to just use xpos (ypos will have the same areas of missing data) |
---|
| 826 | |
---|
| 827 | jimin=1 |
---|
| 828 | jimax=nlons |
---|
| 829 | jjmin=1 |
---|
| 830 | jjmax=nlats |
---|
| 831 | |
---|
| 832 | minlon_xpos: DO ji= 1, nlons |
---|
| 833 | IF (COUNT(ixposi(ji,:) >= 0) > 0) THEN |
---|
| 834 | jimin=ji |
---|
| 835 | EXIT minlon_xpos |
---|
| 836 | ENDIF |
---|
| 837 | END DO minlon_xpos |
---|
| 838 | |
---|
| 839 | maxlon_xpos: DO ji= nlons, 1, -1 |
---|
| 840 | IF (COUNT(ixposi(ji,:) >= 0) > 0) THEN |
---|
| 841 | jimax=ji |
---|
| 842 | EXIT maxlon_xpos |
---|
| 843 | ENDIF |
---|
| 844 | END DO maxlon_xpos |
---|
| 845 | |
---|
| 846 | minlat_xpos: DO jj= 1, nlats |
---|
| 847 | IF (COUNT(ixposi(:,jj) >= 0) > 0) THEN |
---|
| 848 | jjmin=jj |
---|
| 849 | EXIT minlat_xpos |
---|
| 850 | ENDIF |
---|
| 851 | END DO minlat_xpos |
---|
| 852 | |
---|
| 853 | maxlat_xpos: DO jj= nlats, 1, -1 |
---|
| 854 | IF (COUNT(ixposi(:,jj) >= 0) > 0) THEN |
---|
| 855 | jjmax=jj |
---|
| 856 | EXIT maxlat_xpos |
---|
| 857 | ENDIF |
---|
| 858 | END DO maxlat_xpos |
---|
| 859 | |
---|
| 860 | lonmin = lonsi(jimin,jjmin) |
---|
| 861 | latmin = latsi(jimin,jjmin) |
---|
| 862 | nlons = jimax-jimin+1 |
---|
| 863 | nlats = jjmax-jjmin+1 |
---|
| 864 | |
---|
| 865 | ! construct new arrays |
---|
| 866 | |
---|
| 867 | ALLOCATE( & |
---|
| 868 | & lons(nlons,nlats), & |
---|
| 869 | & lats(nlons,nlats), & |
---|
| 870 | & ixpos(nlons,nlats), & |
---|
| 871 | & iypos(nlons,nlats), & |
---|
[2358] | 872 | & iprocn(nlons,nlats) & |
---|
[2128] | 873 | & ) |
---|
| 874 | |
---|
| 875 | lons(:,:) = lonsi(jimin:jimax,jjmin:jjmax) |
---|
| 876 | lats(:,:) = latsi(jimin:jimax,jjmin:jjmax) |
---|
| 877 | ixpos(:,:) = ixposi(jimin:jimax,jjmin:jjmax) |
---|
| 878 | iypos(:,:) = iyposi(jimin:jimax,jjmin:jjmax) |
---|
[2358] | 879 | iprocn(:,:) = iproci(jimin:jimax,jjmin:jjmax) |
---|
[2128] | 880 | |
---|
| 881 | DEALLOCATE(lonsi,latsi,ixposi,iyposi,iproci) |
---|
| 882 | |
---|
| 883 | ! calculate (estimate) maxxdiff, maxydiff |
---|
| 884 | ! this is to help define the search area for obs_grid_search_lookup |
---|
| 885 | |
---|
| 886 | maxxdiff = 1 |
---|
| 887 | maxydiff = 1 |
---|
| 888 | |
---|
| 889 | tmpx1 = 0 |
---|
| 890 | tmpx2 = 0 |
---|
| 891 | tmpy1 = 0 |
---|
| 892 | tmpy2 = 0 |
---|
| 893 | |
---|
| 894 | numx1 = 0 |
---|
| 895 | numx2 = 0 |
---|
| 896 | numy1 = 0 |
---|
| 897 | numy2 = 0 |
---|
| 898 | |
---|
| 899 | ! calculate the mean absolute xdiff and ydiff |
---|
| 900 | ! also calculate a histogram |
---|
| 901 | ! note the reason why looking for xdiff and ydiff in both directions |
---|
| 902 | ! is to allow for rotated grids |
---|
| 903 | |
---|
| 904 | DO ji = 1, nlons-1 |
---|
| 905 | DO jj = 1, nlats-1 |
---|
| 906 | IF ( ixpos(ji,jj) > 0 .AND. iypos(ji,jj) > 0 ) THEN |
---|
| 907 | IF ( ixpos(ji+1,jj) > 0 ) THEN |
---|
| 908 | df = ABS( ixpos(ji+1,jj) - ixpos(ji,jj) ) |
---|
| 909 | tmpx1 = tmpx1+df |
---|
| 910 | numx1 = numx1+1 |
---|
| 911 | IF ( df < histsize ) histx1(df+1) = histx1(df+1) + 1 |
---|
| 912 | ENDIF |
---|
| 913 | IF ( ixpos(ji,jj+1) > 0 ) THEN |
---|
| 914 | df = ABS( ixpos(ji,jj+1) - ixpos(ji,jj) ) |
---|
| 915 | tmpx2 = tmpx2 + df |
---|
| 916 | numx2 = numx2 + 1 |
---|
| 917 | IF ( df < histsize ) histx2(df+1) = histx2(df+1) + 1 |
---|
| 918 | ENDIF |
---|
| 919 | IF (iypos(ji+1,jj) > 0) THEN |
---|
| 920 | df = ABS( iypos(ji+1,jj) - iypos(ji,jj) ) |
---|
| 921 | tmpy1 = tmpy1 + df |
---|
| 922 | numy1 = numy1 + 1 |
---|
| 923 | IF ( df < histsize ) histy1(df+1) = histy1(df+1) + 1 |
---|
| 924 | ENDIF |
---|
| 925 | IF ( iypos(ji,jj+1) > 0 ) THEN |
---|
| 926 | df = ABS( iypos(ji,jj+1) - iypos(ji,jj) ) |
---|
| 927 | tmpy2 = tmpy2 + df |
---|
| 928 | numy2 = numy2 + 1 |
---|
| 929 | IF ( df < histsize ) histy2(df+1) = histy2(df+1) + 1 |
---|
| 930 | ENDIF |
---|
| 931 | ENDIF |
---|
| 932 | END DO |
---|
| 933 | END DO |
---|
| 934 | |
---|
| 935 | IF (lwp) THEN |
---|
| 936 | WRITE(numout,*) 'histograms' |
---|
| 937 | WRITE(numout,*) '0 1 2 3 4 5 6 7 8 9 10 ...' |
---|
| 938 | WRITE(numout,*) 'histx1' |
---|
| 939 | WRITE(numout,*) histx1 |
---|
| 940 | WRITE(numout,*) 'histx2' |
---|
| 941 | WRITE(numout,*) histx2 |
---|
| 942 | WRITE(numout,*) 'histy1' |
---|
| 943 | WRITE(numout,*) histy1 |
---|
| 944 | WRITE(numout,*) 'histy2' |
---|
| 945 | WRITE(numout,*) histy2 |
---|
| 946 | ENDIF |
---|
| 947 | |
---|
| 948 | meanxdiff1 = tmpx1 / numx1 |
---|
| 949 | meanydiff1 = tmpy1 / numy1 |
---|
| 950 | meanxdiff2 = tmpx2 / numx2 |
---|
| 951 | meanydiff2 = tmpy2 / numy2 |
---|
| 952 | |
---|
| 953 | meanxdiff = MAXVAL((/ meanxdiff1, meanxdiff2 /)) |
---|
| 954 | meanydiff = MAXVAL((/ meanydiff1, meanydiff2 /)) |
---|
| 955 | |
---|
| 956 | IF (lwp) THEN |
---|
| 957 | WRITE(numout,*) tmpx1, tmpx2, tmpy1, tmpy2 |
---|
| 958 | WRITE(numout,*) numx1, numx2, numy1, numy2 |
---|
| 959 | WRITE(numout,*) 'meanxdiff: ',meanxdiff, meanxdiff1, meanxdiff2 |
---|
| 960 | WRITE(numout,*) 'meanydiff: ',meanydiff, meanydiff1, meanydiff2 |
---|
| 961 | ENDIF |
---|
| 962 | |
---|
| 963 | tmpx1 = 0 |
---|
| 964 | tmpx2 = 0 |
---|
| 965 | tmpy1 = 0 |
---|
| 966 | tmpy2 = 0 |
---|
| 967 | |
---|
| 968 | numx1 = 0 |
---|
| 969 | numx2 = 0 |
---|
| 970 | numy1 = 0 |
---|
| 971 | numy2 = 0 |
---|
| 972 | |
---|
| 973 | histx1(:) = 0 |
---|
| 974 | histx2(:) = 0 |
---|
| 975 | histy1(:) = 0 |
---|
| 976 | histy2(:) = 0 |
---|
| 977 | |
---|
| 978 | limxdiff = meanxdiff * 4! limit the difference to avoid picking up wraparound |
---|
| 979 | limydiff = meanydiff * 4 |
---|
| 980 | |
---|
| 981 | DO ji = 1, nlons-1 |
---|
| 982 | DO jj = 1, nlats-1 |
---|
| 983 | IF ( ixpos(ji,jj) > 0 .AND. iypos(ji,jj) > 0 ) THEN |
---|
| 984 | |
---|
| 985 | IF ( ixpos(ji+1,jj) > 0 ) THEN |
---|
| 986 | df = ABS( ixpos(ji+1,jj)-ixpos(ji,jj) ) |
---|
| 987 | tmpx1 = df |
---|
| 988 | IF ( df < limxdiff ) numx1 = numx1+1 |
---|
| 989 | IF ( df < histsize ) histx1(df+1) = histx1(df+1) + 1 |
---|
| 990 | ENDIF |
---|
| 991 | IF ( ixpos(ji,jj+1) > 0 ) THEN |
---|
| 992 | df = ABS( ixpos(ji,jj+1) - ixpos(ji,jj) ) |
---|
| 993 | tmpx2 = df |
---|
| 994 | IF ( df < limxdiff ) numx2 = numx2 + 1 |
---|
| 995 | IF ( df < histsize ) histx2(df+1) = histx2(df+1) + 1 |
---|
| 996 | ENDIF |
---|
| 997 | IF (iypos(ji+1,jj) > 0) THEN |
---|
| 998 | df = ABS( iypos(ji+1,jj) - iypos(ji,jj) ) |
---|
| 999 | tmpy1 = df |
---|
[4990] | 1000 | IF ( df < limydiff ) numy1 = numy1 + 1 |
---|
[2128] | 1001 | IF ( df < histsize ) histy1(df+1) = histy1(df+1) + 1 |
---|
| 1002 | ENDIF |
---|
| 1003 | IF (iypos(ji,jj+1) > 0) THEN |
---|
| 1004 | df = ABS( iypos(ji,jj+1) - iypos(ji,jj) ) |
---|
| 1005 | tmpy2 = df |
---|
[4990] | 1006 | IF ( df < limydiff ) numy2 = numy2+1 |
---|
[2128] | 1007 | IF ( df < histsize ) histy2(df+1) = histy2(df+1)+1 |
---|
| 1008 | ENDIF |
---|
| 1009 | |
---|
| 1010 | IF ( maxxdiff < tmpx1 .AND. tmpx1 < limxdiff ) & |
---|
| 1011 | & maxxdiff = tmpx1 |
---|
| 1012 | IF ( maxxdiff < tmpx2 .AND. tmpx2 < limxdiff ) & |
---|
| 1013 | & maxxdiff = tmpx2 |
---|
| 1014 | IF ( maxydiff < tmpy1 .AND. tmpy1 < limydiff ) & |
---|
| 1015 | & maxydiff = tmpy1 |
---|
| 1016 | IF ( maxydiff < tmpy2 .AND. tmpy2 < limydiff ) & |
---|
| 1017 | & maxydiff = tmpy2 |
---|
| 1018 | |
---|
| 1019 | ENDIF |
---|
| 1020 | END DO |
---|
| 1021 | END DO |
---|
| 1022 | |
---|
| 1023 | ! cumulative histograms |
---|
| 1024 | |
---|
| 1025 | DO ji = 1, histsize - 1 |
---|
| 1026 | histx1(ji+1) = histx1(ji+1) + histx1(ji) |
---|
| 1027 | histx2(ji+1) = histx2(ji+1) + histx2(ji) |
---|
| 1028 | histy1(ji+1) = histy1(ji+1) + histy1(ji) |
---|
| 1029 | histy2(ji+1) = histy2(ji+1) + histy2(ji) |
---|
| 1030 | END DO |
---|
| 1031 | |
---|
| 1032 | fhistx1(:) = histx1(:) * 1.0 / numx1 |
---|
| 1033 | fhistx2(:) = histx2(:) * 1.0 / numx2 |
---|
| 1034 | fhisty1(:) = histy1(:) * 1.0 / numy1 |
---|
| 1035 | fhisty2(:) = histy2(:) * 1.0 / numy2 |
---|
| 1036 | |
---|
| 1037 | ! output new histograms |
---|
| 1038 | |
---|
| 1039 | IF (lwp) THEN |
---|
| 1040 | WRITE(numout,*) 'cumulative histograms' |
---|
| 1041 | WRITE(numout,*) '0 1 2 3 4 5 6 7 8 9 10 ...' |
---|
| 1042 | WRITE(numout,*) 'fhistx1' |
---|
| 1043 | WRITE(numout,*) fhistx1 |
---|
| 1044 | WRITE(numout,*) 'fhistx2' |
---|
| 1045 | WRITE(numout,*) fhistx2 |
---|
| 1046 | WRITE(numout,*) 'fhisty1' |
---|
| 1047 | WRITE(numout,*) fhisty1 |
---|
| 1048 | WRITE(numout,*) 'fhisty2' |
---|
| 1049 | WRITE(numout,*) fhisty2 |
---|
| 1050 | ENDIF |
---|
| 1051 | |
---|
| 1052 | ! calculate maxxdiff and maxydiff based on cumulative histograms |
---|
| 1053 | ! where > 0.999 of points are |
---|
| 1054 | |
---|
| 1055 | ! maxval just converts 1x1 vector return from maxloc to a scalar |
---|
| 1056 | |
---|
| 1057 | histtol = 0.999 |
---|
| 1058 | tmpx1 = MAXVAL( MAXLOC( fhistx1(:), mask = ( fhistx1(:) <= histtol ) ) ) |
---|
| 1059 | tmpx2 = MAXVAL( MAXLOC( fhistx2(:), mask = ( fhistx2(:) <= histtol ) ) ) |
---|
| 1060 | tmpy1 = MAXVAL( MAXLOC( fhisty1(:), mask = ( fhisty1(:) <= histtol ) ) ) |
---|
| 1061 | tmpy2 = MAXVAL( MAXLOC( fhisty2(:), mask = ( fhisty2(:) <= histtol ) ) ) |
---|
| 1062 | |
---|
| 1063 | maxxdiff = MAXVAL( (/ tmpx1, tmpx2 /) ) + 1 |
---|
| 1064 | maxydiff = MAXVAL( (/ tmpy1, tmpy2 /) ) + 1 |
---|
| 1065 | |
---|
| 1066 | ! Write out data |
---|
| 1067 | |
---|
| 1068 | IF ( ( .NOT. ln_grid_global ) .OR. & |
---|
| 1069 | & ( ( ln_grid_global ) .AND. ( nproc==0 ) ) ) THEN |
---|
| 1070 | |
---|
| 1071 | CALL chkerr( nf90_create (TRIM(cfname), nf90_clobber, idfile), & |
---|
| 1072 | & cpname, __LINE__ ) |
---|
| 1073 | CALL chkerr( nf90_put_att( idfile, nf90_global, 'title', & |
---|
| 1074 | & 'Mapping file from lon/lat to model grid point' ),& |
---|
| 1075 | & cpname,__LINE__ ) |
---|
| 1076 | CALL chkerr( nf90_put_att( idfile, nf90_global, 'maxxdiff', & |
---|
| 1077 | & maxxdiff ), & |
---|
| 1078 | & cpname,__LINE__ ) |
---|
| 1079 | CALL chkerr( nf90_put_att( idfile, nf90_global, 'maxydiff', & |
---|
| 1080 | & maxydiff ), & |
---|
| 1081 | & cpname,__LINE__ ) |
---|
| 1082 | CALL chkerr( nf90_put_att( idfile, nf90_global, 'dlon', dlon ),& |
---|
| 1083 | & cpname,__LINE__ ) |
---|
| 1084 | CALL chkerr( nf90_put_att( idfile, nf90_global, 'dlat', dlat ),& |
---|
| 1085 | & cpname,__LINE__ ) |
---|
| 1086 | CALL chkerr( nf90_put_att( idfile, nf90_global, 'lonmin', & |
---|
| 1087 | & lonmin ), & |
---|
| 1088 | & cpname,__LINE__ ) |
---|
| 1089 | CALL chkerr( nf90_put_att( idfile, nf90_global, 'latmin', & |
---|
| 1090 | & latmin ), & |
---|
| 1091 | & cpname,__LINE__ ) |
---|
| 1092 | |
---|
| 1093 | CALL chkerr( nf90_def_dim(idfile, 'nx' , nlons, idnx), & |
---|
| 1094 | & cpname,__LINE__ ) |
---|
| 1095 | CALL chkerr( nf90_def_dim(idfile, 'ny' , nlats, idny), & |
---|
| 1096 | & cpname,__LINE__ ) |
---|
| 1097 | |
---|
| 1098 | incdim(1) = idnx |
---|
| 1099 | incdim(2) = idny |
---|
| 1100 | |
---|
| 1101 | CALL chkerr( nf90_def_var( idfile, 'LON', nf90_float, incdim, & |
---|
| 1102 | & idlon ), & |
---|
| 1103 | & cpname, __LINE__ ) |
---|
| 1104 | CALL chkerr( nf90_put_att( idfile, idlon, 'long_name', & |
---|
| 1105 | & 'longitude' ), & |
---|
| 1106 | & cpname, __LINE__ ) |
---|
| 1107 | |
---|
| 1108 | CALL chkerr( nf90_def_var( idfile, 'LAT', nf90_float, incdim, & |
---|
| 1109 | & idlat ), & |
---|
| 1110 | & cpname, __LINE__ ) |
---|
| 1111 | CALL chkerr( nf90_put_att( idfile, idlat, 'long_name', & |
---|
| 1112 | & 'latitude' ), & |
---|
| 1113 | & cpname, __LINE__ ) |
---|
| 1114 | |
---|
| 1115 | CALL chkerr( nf90_def_var( idfile, 'XPOS', nf90_int, incdim, & |
---|
| 1116 | & idxpos ), & |
---|
| 1117 | & cpname, __LINE__ ) |
---|
| 1118 | CALL chkerr( nf90_put_att( idfile, idxpos, 'long_name', & |
---|
| 1119 | & 'x position' ), & |
---|
| 1120 | & cpname, __LINE__ ) |
---|
| 1121 | CALL chkerr( nf90_put_att( idfile, idxpos, '_FillValue', -1 ), & |
---|
| 1122 | & cpname, __LINE__ ) |
---|
| 1123 | |
---|
| 1124 | CALL chkerr( nf90_def_var( idfile, 'YPOS', nf90_int, incdim, & |
---|
| 1125 | & idypos ), & |
---|
| 1126 | & cpname, __LINE__ ) |
---|
| 1127 | CALL chkerr( nf90_put_att( idfile, idypos, 'long_name', & |
---|
| 1128 | & 'y position' ), & |
---|
| 1129 | & cpname, __LINE__ ) |
---|
| 1130 | CALL chkerr( nf90_put_att( idfile, idypos, '_FillValue', -1 ), & |
---|
| 1131 | & cpname, __LINE__ ) |
---|
| 1132 | |
---|
| 1133 | CALL chkerr( nf90_enddef( idfile ), cpname, __LINE__ ) |
---|
| 1134 | |
---|
| 1135 | CALL chkerr( nf90_put_var( idfile, idlon, lons), & |
---|
| 1136 | & cpname, __LINE__ ) |
---|
| 1137 | CALL chkerr( nf90_put_var( idfile, idlat, lats), & |
---|
| 1138 | & cpname, __LINE__ ) |
---|
| 1139 | CALL chkerr( nf90_put_var( idfile, idxpos, ixpos), & |
---|
| 1140 | & cpname, __LINE__ ) |
---|
| 1141 | CALL chkerr( nf90_put_var( idfile, idypos, iypos), & |
---|
| 1142 | & cpname, __LINE__ ) |
---|
| 1143 | |
---|
| 1144 | CALL chkerr( nf90_close( idfile ), cpname, __LINE__ ) |
---|
| 1145 | |
---|
| 1146 | ! should also output max i, max j spacing for use in |
---|
| 1147 | ! obs_grid_search_lookup |
---|
| 1148 | |
---|
| 1149 | ENDIF |
---|
| 1150 | |
---|
| 1151 | ENDIF |
---|
| 1152 | |
---|
| 1153 | ENDIF |
---|
| 1154 | |
---|
| 1155 | END SUBROUTINE obs_grid_setup |
---|
| 1156 | |
---|
| 1157 | SUBROUTINE obs_grid_deallocate( ) |
---|
| 1158 | !!---------------------------------------------------------------------- |
---|
| 1159 | !! *** ROUTINE obs_grid_setup *** |
---|
| 1160 | !! |
---|
| 1161 | !! ** Purpose : Deallocate arrays setup by obs_grid_setup |
---|
| 1162 | !! |
---|
| 1163 | !! History : |
---|
| 1164 | !! ! 2007-12 (D. Lea) new routine |
---|
| 1165 | !!----------------------------------------------------------------------- |
---|
| 1166 | |
---|
| 1167 | IF (ln_grid_search_lookup) THEN |
---|
[2358] | 1168 | DEALLOCATE( lons, lats, ixpos, iypos, iprocn ) |
---|
[2128] | 1169 | ENDIF |
---|
| 1170 | |
---|
| 1171 | END SUBROUTINE obs_grid_deallocate |
---|
| 1172 | |
---|
| 1173 | #include "obs_level_search.h90" |
---|
| 1174 | |
---|
| 1175 | #include "linquad.h90" |
---|
| 1176 | |
---|
| 1177 | #include "maxdist.h90" |
---|
| 1178 | |
---|
| 1179 | #include "find_obs_proc.h90" |
---|
| 1180 | |
---|
| 1181 | END MODULE obs_grid |
---|
| 1182 | |
---|