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