[8] | 1 | |
---|
| 2 | !! This module define variables for the grid to gathered points. |
---|
| 3 | !! |
---|
| 4 | !! @call sechiba_main |
---|
[12] | 5 | !! @Version : $Revision$, $Date$ |
---|
[8] | 6 | !! |
---|
[12] | 7 | !< $HeadURL$ |
---|
| 8 | !< $Date$ |
---|
| 9 | !< $Revision$ |
---|
[8] | 10 | !! |
---|
| 11 | !! @author Marie-Alice Foujols, Jan Polcher and Martial Mancip |
---|
| 12 | !! |
---|
[3447] | 13 | !! This module archives and makes available for all ORCHIDEE routine the information on the grid |
---|
| 14 | !! being used. 3 types of grids are foreseen : |
---|
| 15 | !! - Regular longitude latitude grid : This is the default and mostly used for global applications. |
---|
| 16 | !! - Regular X/Y grid : this is a typical grid for regional models and requires a projection method |
---|
| 17 | !! to go from X/y to lon/lat. |
---|
[5559] | 18 | !! - unstructures grid : This is a general grid where each cell is a polygone. It is used for DYNAMICO. |
---|
[3447] | 19 | !! |
---|
| 20 | !! The subroutines have the following role : |
---|
| 21 | !! grid_init : this routine will provide the dimensions needed to allocate the memory and the |
---|
| 22 | !! characteristics of the grid. |
---|
| 23 | !! |
---|
| 24 | !! grid_stuff : This subroutine provides the grid details for all land points. Obviously depending |
---|
| 25 | !! on the grid type different level of information need to be provided. |
---|
| 26 | !! |
---|
[8] | 27 | MODULE grid |
---|
| 28 | |
---|
[4179] | 29 | USE grid_var |
---|
[8] | 30 | USE defprec |
---|
[5364] | 31 | USE mod_orchidee_para_var |
---|
| 32 | USE mod_orchidee_transfert_para |
---|
[3447] | 33 | USE haversine |
---|
| 34 | USE module_llxy |
---|
| 35 | |
---|
| 36 | USE ioipsl |
---|
| 37 | USE netcdf |
---|
| 38 | |
---|
[8] | 39 | IMPLICIT NONE |
---|
[5559] | 40 | |
---|
| 41 | PUBLIC grid_init, grid_init_unstructured, grid_set_glo, grid_allocate_glo, grid_stuff, grid_tolola, grid_initproj |
---|
| 42 | |
---|
[8] | 43 | ! |
---|
[3447] | 44 | !================================================================================= |
---|
| 45 | ! |
---|
| 46 | ! Horizontal grid information |
---|
| 47 | ! |
---|
| 48 | !================================================================================= |
---|
[4179] | 49 | |
---|
[3447] | 50 | ! Global map or not. |
---|
| 51 | ! There is little chance that if iim <=2 and jjm <= 2 that we have global grid. |
---|
| 52 | ! Furthermore using the second line allows to avoid pole problems for global grids |
---|
| 53 | LOGICAL, SAVE :: global = .TRUE. |
---|
| 54 | !$OMP THREADPRIVATE(global) |
---|
[4179] | 55 | |
---|
[8] | 56 | ! PARAMETERS |
---|
[5364] | 57 | REAL(r_std), PARAMETER :: default_resolution = 250000. !! default resolution (m) |
---|
| 58 | |
---|
[8] | 59 | ! VARIABLES |
---|
| 60 | ! |
---|
| 61 | !- |
---|
| 62 | !- Variable to help describe the grid |
---|
| 63 | !- once the points are gathered. |
---|
| 64 | !- |
---|
[5364] | 65 | REAL(r_std), SAVE :: limit_west, limit_east !! Limits of the domain |
---|
| 66 | REAL(r_std), SAVE :: limit_north, limit_south !! Limits of the domain |
---|
[1078] | 67 | !$OMP THREADPRIVATE(limit_west, limit_east, limit_north, limit_south) |
---|
[8] | 68 | !- |
---|
| 69 | !! Geographical coordinates |
---|
| 70 | REAL(r_std), ALLOCATABLE, DIMENSION (:,:), SAVE :: lalo |
---|
[1078] | 71 | !$OMP THREADPRIVATE(lalo) |
---|
[8] | 72 | !! index of land points |
---|
[3447] | 73 | INTEGER, ALLOCATABLE, DIMENSION (:), SAVE :: ilandindex,jlandindex |
---|
[1078] | 74 | !$OMP THREADPRIVATE(ilandindex, jlandindex) |
---|
[8] | 75 | !- |
---|
| 76 | !! Fraction of continents. |
---|
| 77 | REAL(r_std), ALLOCATABLE, DIMENSION (:), SAVE :: contfrac |
---|
[1078] | 78 | !$OMP THREADPRIVATE(contfrac) |
---|
[8] | 79 | ! |
---|
[3447] | 80 | ! indices of the NbNeighb neighbours of each grid point |
---|
| 81 | ! (1=Northern most vertex and then in clockwise order) |
---|
| 82 | ! Zero or negative index means that this neighbour is not a land point |
---|
| 83 | INTEGER(i_std), ALLOCATABLE, DIMENSION (:,:), SAVE :: neighbours |
---|
[1078] | 84 | !$OMP THREADPRIVATE(neighbours) |
---|
[8] | 85 | ! |
---|
[3447] | 86 | ! Heading of the direction out of the grid box either through the vertex |
---|
| 87 | ! of the mid-segment of the polygon. |
---|
| 88 | ! |
---|
| 89 | REAL(r_std), ALLOCATABLE, DIMENSION(:,:), SAVE :: headings |
---|
| 90 | !$OMP THREADPRIVATE(headings) |
---|
| 91 | ! |
---|
| 92 | ! Length of segments of the polygon. |
---|
| 93 | ! |
---|
| 94 | REAL(r_std), ALLOCATABLE, DIMENSION(:,:), SAVE :: seglength |
---|
| 95 | !$OMP THREADPRIVATE(seglength) |
---|
| 96 | ! |
---|
| 97 | ! Area of the grid box |
---|
| 98 | ! |
---|
| 99 | REAL(r_std), ALLOCATABLE, DIMENSION(:), SAVE :: area |
---|
| 100 | !$OMP THREADPRIVATE(area) |
---|
| 101 | ! |
---|
| 102 | ! Coordinats of the vertices |
---|
| 103 | ! |
---|
| 104 | REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:), SAVE :: corners |
---|
| 105 | !$OMP THREADPRIVATE(corners) |
---|
| 106 | ! |
---|
| 107 | ! Resolution remains a temporary variable until the merge of the |
---|
| 108 | ! re-interfacing of the interpolation by Lluis. One this is done |
---|
| 109 | ! Resolution will be replaced in the model either by area or seglength. |
---|
| 110 | ! |
---|
[8] | 111 | REAL(r_std), ALLOCATABLE, DIMENSION (:,:), SAVE :: resolution |
---|
[1078] | 112 | !$OMP THREADPRIVATE(resolution) |
---|
[8] | 113 | ! |
---|
| 114 | ! |
---|
[3447] | 115 | ! |
---|
| 116 | INTEGER(i_std), PARAMETER :: MAX_DOMAINS=1 |
---|
[8] | 117 | ! |
---|
[3447] | 118 | type (proj_info), SAVE, dimension(1:MAX_DOMAINS) :: proj_stack |
---|
| 119 | ! |
---|
| 120 | real(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: dxwrf, dywrf |
---|
| 121 | ! |
---|
| 122 | ! |
---|
| 123 | INTERFACE grid_tolola |
---|
| 124 | MODULE PROCEDURE grid_tolola_scal, grid_tolola_1d, grid_tolola_2d |
---|
| 125 | END INTERFACE grid_tolola |
---|
| 126 | |
---|
| 127 | INTERFACE grid_toij |
---|
| 128 | MODULE PROCEDURE grid_toij_scal, grid_toij_1d, grid_toij_2d |
---|
| 129 | END INTERFACE grid_toij |
---|
| 130 | ! |
---|
[8] | 131 | CONTAINS |
---|
[5559] | 132 | |
---|
[3447] | 133 | !! ============================================================================================================================= |
---|
| 134 | !! SUBROUTINE: grid_init |
---|
| 135 | !! |
---|
| 136 | !>\BRIEF Initialization of grid description distributed by this module to the rest of the model. |
---|
| 137 | !! |
---|
| 138 | !! DESCRIPTION: Routine which provides the dimension of the grid (number of land points) as well as the |
---|
| 139 | !! grid characteristics (type and name) so that the memory can be allocated. |
---|
| 140 | !! |
---|
[5364] | 141 | !! This subroutine is called by intersurf_main_gathered, readdim2, orchideedriver or any driver of the model. |
---|
[3447] | 142 | !! |
---|
| 143 | !! \n |
---|
| 144 | !_ ============================================================================================================================== |
---|
| 145 | !! |
---|
| 146 | SUBROUTINE grid_init ( npts, nbseg, gtype, gname, isglobal ) |
---|
[8] | 147 | ! |
---|
| 148 | ! 0 interface |
---|
| 149 | ! |
---|
| 150 | IMPLICIT NONE |
---|
| 151 | ! |
---|
| 152 | ! 0.1 input ! |
---|
[3447] | 153 | ! |
---|
[8] | 154 | ! Domain size |
---|
[3447] | 155 | INTEGER(i_std), INTENT(in) :: npts !! Number of local continental points |
---|
| 156 | INTEGER(i_std), INTENT(in) :: nbseg !! number of segments of the polygone of the mesh |
---|
[5559] | 157 | INTEGER(i_std), INTENT(in) :: gtype !! Type of grid: regular_lonlat/regular_xy/unstructured |
---|
[3447] | 158 | CHARACTER(LEN=*), INTENT(in) :: gname !! Name of the grid |
---|
| 159 | LOGICAL, OPTIONAL :: isglobal |
---|
[8] | 160 | ! |
---|
[3447] | 161 | ! Verify the information passed and save it in the global variables of the model. |
---|
| 162 | ! |
---|
[5559] | 163 | IF ( gtype == regular_lonlat ) THEN |
---|
| 164 | grid_type=regular_lonlat |
---|
[3447] | 165 | |
---|
[5364] | 166 | IF ( nbseg /= 4 ) THEN |
---|
[3447] | 167 | CALL ipslerr(3, "grid_init", "This regular Lon/lat grid should have 4 segments", & |
---|
| 168 | & "per horizontal grid box","") |
---|
| 169 | ELSE |
---|
| 170 | NbSegments=4 |
---|
| 171 | ENDIF |
---|
| 172 | IF ( PRESENT(isglobal) ) THEN |
---|
| 173 | global = isglobal |
---|
| 174 | ELSE |
---|
| 175 | global = .TRUE. |
---|
| 176 | ENDIF |
---|
[5364] | 177 | |
---|
[5559] | 178 | ELSE IF ( gtype == regular_xy ) THEN |
---|
| 179 | grid_type=regular_xy |
---|
[3447] | 180 | IF ( nbseg /= 4 ) THEN |
---|
| 181 | CALL ipslerr(3, "grid_init", "This regular X/Y grid should have 4 segments", & |
---|
| 182 | & "per horizontal grid box","") |
---|
| 183 | ELSE |
---|
| 184 | NbSegments=4 |
---|
| 185 | ENDIF |
---|
| 186 | IF ( PRESENT(isglobal) ) THEN |
---|
| 187 | global = isglobal |
---|
| 188 | ELSE |
---|
| 189 | global = .FALSE. |
---|
| 190 | ENDIF |
---|
[5559] | 191 | ELSE IF ( gtype == unstructured ) THEN |
---|
| 192 | grid_type=unstructured |
---|
[3447] | 193 | NbSegments=nbseg |
---|
| 194 | IF ( PRESENT(isglobal) ) THEN |
---|
| 195 | global = isglobal |
---|
| 196 | ELSE |
---|
| 197 | global = .TRUE. |
---|
| 198 | ENDIF |
---|
| 199 | ELSE |
---|
[5559] | 200 | WRITE(numout,*) "Unrecognized grid type: gtype=",gtype |
---|
[3447] | 201 | CALL ipslerr(3, "grid_init", "unrecognized grid type.",& |
---|
[5559] | 202 | & "It has to be either regular_lonlat, regular_xy or unstructured","") |
---|
[3447] | 203 | ENDIF |
---|
| 204 | ! |
---|
[8] | 205 | ! Create the internal coordinate table |
---|
| 206 | ! |
---|
| 207 | IF ( (.NOT.ALLOCATED(lalo))) THEN |
---|
| 208 | ALLOCATE(lalo(npts,2)) |
---|
| 209 | lalo(:,:) = val_exp |
---|
| 210 | ENDIF |
---|
| 211 | !- |
---|
| 212 | !- Store variable to help describe the grid |
---|
| 213 | !- once the points are gathered. |
---|
| 214 | !- |
---|
[3447] | 215 | NbNeighb=2*NbSegments |
---|
[8] | 216 | IF ( (.NOT.ALLOCATED(neighbours))) THEN |
---|
[3447] | 217 | ALLOCATE(neighbours(npts,NbNeighb)) |
---|
[8] | 218 | neighbours(:,:) = -999999 |
---|
| 219 | ENDIF |
---|
[3447] | 220 | IF ( (.NOT.ALLOCATED(headings))) THEN |
---|
| 221 | ALLOCATE(headings(npts,NbNeighb)) |
---|
| 222 | headings(:,:) = val_exp |
---|
[8] | 223 | ENDIF |
---|
[3447] | 224 | IF ( (.NOT.ALLOCATED(seglength))) THEN |
---|
| 225 | ALLOCATE(seglength(npts,NbSegments)) |
---|
| 226 | seglength(:,:) = val_exp |
---|
| 227 | ENDIF |
---|
| 228 | IF ( (.NOT.ALLOCATED(corners))) THEN |
---|
| 229 | ALLOCATE(corners(npts,NbSegments,2)) |
---|
| 230 | corners(:,:,:) = val_exp |
---|
| 231 | ENDIF |
---|
[8] | 232 | IF ( (.NOT.ALLOCATED(area))) THEN |
---|
| 233 | ALLOCATE(area(npts)) |
---|
| 234 | area(:) = val_exp |
---|
| 235 | ENDIF |
---|
| 236 | ! |
---|
[3447] | 237 | ! TEMPORARY |
---|
| 238 | ! |
---|
| 239 | IF ( (.NOT.ALLOCATED(resolution))) THEN |
---|
| 240 | ALLOCATE(resolution(npts,2)) |
---|
| 241 | resolution(:,:) = val_exp |
---|
| 242 | ENDIF |
---|
| 243 | ! |
---|
[8] | 244 | !- Store the fraction of the continents only once so that the user |
---|
| 245 | !- does not change them afterwards. |
---|
| 246 | ! |
---|
| 247 | IF ( (.NOT.ALLOCATED(contfrac))) THEN |
---|
| 248 | ALLOCATE(contfrac(npts)) |
---|
| 249 | contfrac(:) = val_exp |
---|
| 250 | ENDIF |
---|
| 251 | ! |
---|
[3447] | 252 | ! Allocation of index coordinates ... |
---|
| 253 | ! JP : these are global fields and should perhaps be allocated somewhere else. |
---|
[8] | 254 | IF (.NOT. ALLOCATED(ilandindex)) THEN |
---|
[3447] | 255 | ALLOCATE(ilandindex(nbp_glo),jlandindex(nbp_glo)) |
---|
[8] | 256 | ilandindex(:) = -10000000 |
---|
| 257 | jlandindex(:) = -10000000 |
---|
| 258 | ENDIF |
---|
| 259 | ! |
---|
[3447] | 260 | END SUBROUTINE grid_init |
---|
[5364] | 261 | |
---|
| 262 | |
---|
| 263 | !! ============================================================================================================================= |
---|
| 264 | !! SUBROUTINE: grid_init_unstructured |
---|
[3447] | 265 | !! |
---|
[5364] | 266 | !>\BRIEF Initialization of grid description for unstructured grid. |
---|
[3447] | 267 | !! |
---|
[5364] | 268 | !! DESCRIPTION: |
---|
| 269 | !! This subroutine is called by intersurf_main_gathered instead of grid_init. grid_init is called in here. |
---|
| 270 | !! |
---|
| 271 | !! \n |
---|
| 272 | !_ ============================================================================================================================== |
---|
[5559] | 273 | SUBROUTINE grid_init_unstructured ( npts, ncells_in, longitude_in, latitude_in, & |
---|
| 274 | bounds_lon_in, bounds_lat_in, area_in, ind_cell_glo_in) |
---|
[5364] | 275 | |
---|
| 276 | IMPLICIT NONE |
---|
| 277 | |
---|
| 278 | INTEGER(i_std), INTENT(IN) :: npts |
---|
[5559] | 279 | INTEGER(i_std), INTENT(IN) :: ncells_in |
---|
| 280 | REAL(r_std), INTENT(IN) :: longitude_in(:) |
---|
| 281 | REAL(r_std), INTENT(IN) :: latitude_in(:) |
---|
| 282 | REAL(r_std), INTENT(IN) :: bounds_lon_in(:,:) |
---|
| 283 | REAL(r_std), INTENT(IN) :: bounds_lat_in(:,:) |
---|
| 284 | REAL(r_std), INTENT(IN) :: area_in(:) |
---|
| 285 | INTEGER(i_std), INTENT(IN) :: ind_cell_glo_in(:) |
---|
[5364] | 286 | |
---|
| 287 | |
---|
[5559] | 288 | CALL grid_init(npts, 6, unstructured, "2DGrid") |
---|
[5364] | 289 | |
---|
| 290 | grid_type=unstructured |
---|
[5559] | 291 | ncells=ncells_in |
---|
[5364] | 292 | nvertex=6 ! for now for dynamico |
---|
| 293 | ALLOCATE(longitude(ncells)) |
---|
| 294 | ALLOCATE(latitude(ncells)) |
---|
| 295 | ALLOCATE(bounds_lon(ncells,nvertex)) |
---|
| 296 | ALLOCATE(bounds_lat(ncells,nvertex)) |
---|
| 297 | ALLOCATE(ind_cell_glo(ncells)) |
---|
| 298 | |
---|
[5559] | 299 | longitude(:) = longitude_in(:) |
---|
| 300 | latitude(:) = latitude_in(:) |
---|
| 301 | bounds_lon(:,:) = bounds_lon_in(:,:) |
---|
| 302 | bounds_lat(:,:) = bounds_lat_in(:,:) |
---|
| 303 | area(:) = area_in(:) |
---|
| 304 | ind_cell_glo(:) = ind_cell_glo_in(:) |
---|
[5364] | 305 | |
---|
| 306 | END SUBROUTINE grid_init_unstructured |
---|
| 307 | |
---|
| 308 | |
---|
[3447] | 309 | !! ============================================================================================================================= |
---|
| 310 | !! FUNCTION grid_set |
---|
| 311 | !! |
---|
| 312 | !>\BRIEF subroutine to set global grid parameters present on all procs |
---|
| 313 | !! |
---|
| 314 | !! DESCRIPTION: |
---|
| 315 | !! |
---|
| 316 | !! |
---|
| 317 | !! |
---|
| 318 | !! |
---|
| 319 | !! \n |
---|
| 320 | !_ ============================================================================================================================== |
---|
| 321 | !! |
---|
| 322 | SUBROUTINE grid_set_glo(arg_nbp_lon,arg_nbp_lat,arg_nbp_glo) |
---|
| 323 | IMPLICIT NONE |
---|
[8] | 324 | |
---|
[3447] | 325 | INTEGER(i_std), INTENT(IN) :: arg_nbp_lon |
---|
| 326 | INTEGER(i_std), INTENT(IN) :: arg_nbp_lat |
---|
| 327 | INTEGER(i_std), INTENT(IN),OPTIONAL :: arg_nbp_glo |
---|
| 328 | iim_g=arg_nbp_lon |
---|
| 329 | jjm_g=arg_nbp_lat |
---|
| 330 | IF (PRESENT(arg_nbp_glo)) nbp_glo=arg_nbp_glo |
---|
| 331 | END SUBROUTINE grid_set_glo |
---|
| 332 | !! ============================================================================================================================= |
---|
| 333 | !! FUNCTION grid_set/allocate_glo |
---|
| 334 | !! |
---|
| 335 | !>\BRIEF subroutines to allocate variables present on all procs |
---|
| 336 | !! |
---|
| 337 | !! DESCRIPTION: |
---|
| 338 | !! |
---|
| 339 | !! |
---|
| 340 | !! |
---|
| 341 | !! |
---|
| 342 | !! \n |
---|
| 343 | !_ ============================================================================================================================== |
---|
| 344 | !! |
---|
| 345 | SUBROUTINE grid_allocate_glo(nbseg) |
---|
[8] | 346 | ! |
---|
[3447] | 347 | IMPLICIT NONE |
---|
| 348 | ! 0.1 input ! |
---|
| 349 | ! |
---|
| 350 | ! Domain size |
---|
| 351 | INTEGER(i_std), INTENT(in) :: nbseg !! number of segments of the polygone of the mesh |
---|
| 352 | ! |
---|
| 353 | ! In case the allocation of the grid is called before the initialisation, |
---|
| 354 | ! we already set the number of segments. |
---|
| 355 | ! This will be done properly in grid_init. |
---|
| 356 | ! |
---|
| 357 | IF ( NbSegments < 3 ) THEN |
---|
| 358 | NbSegments = nbseg |
---|
| 359 | NbNeighb=2*NbSegments |
---|
| 360 | ENDIF |
---|
[6289] | 361 | |
---|
| 362 | |
---|
| 363 | ! Do following only for the OpenMP master thread |
---|
| 364 | IF (is_omp_root) THEN |
---|
| 365 | ALLOCATE(neighbours_g(nbp_glo,NbNeighb)) |
---|
| 366 | ALLOCATE(headings_g(nbp_glo,NbNeighb)) |
---|
| 367 | ALLOCATE(seglength_g(nbp_glo,NbSegments)) |
---|
| 368 | ALLOCATE(corners_g(nbp_glo,NbSegments,2)) |
---|
| 369 | ALLOCATE(area_g(nbp_glo)) |
---|
| 370 | ! |
---|
| 371 | ! TEMPORARY |
---|
| 372 | ! |
---|
| 373 | ALLOCATE(resolution_g(nbp_glo,2)) |
---|
| 374 | ! |
---|
| 375 | ! Allocate other variables |
---|
| 376 | ! |
---|
| 377 | ALLOCATE(lalo_g(nbp_glo,2), contfrac_g(nbp_glo),index_g(nbp_glo)) |
---|
| 378 | ALLOCATE(lon_g(iim_g, jjm_g), lat_g(iim_g, jjm_g)) |
---|
| 379 | END IF |
---|
| 380 | |
---|
[3447] | 381 | END SUBROUTINE grid_allocate_glo |
---|
| 382 | !! |
---|
| 383 | !! ============================================================================================================================= |
---|
| 384 | !! SUBROUTINE: grid_stuff |
---|
| 385 | !! |
---|
| 386 | !>\BRIEF transfers the global horizontal grid information to ORCHIDEE in the case of grid regular in Longitude |
---|
| 387 | !! and Latitude. |
---|
| 388 | !! |
---|
| 389 | !! DESCRIPTION: |
---|
| 390 | !! |
---|
| 391 | !! |
---|
| 392 | !! This subroutine is called by intersurf_main_2d or any driver of the model. |
---|
| 393 | !! |
---|
| 394 | !! \n |
---|
| 395 | !_ ============================================================================================================================== |
---|
| 396 | !! |
---|
| 397 | SUBROUTINE grid_stuff (npts_glo, iim, jjm, grid_lon, grid_lat, kindex, contfrac_tmp) |
---|
| 398 | ! |
---|
[8] | 399 | ! 0 interface |
---|
| 400 | ! |
---|
| 401 | IMPLICIT NONE |
---|
| 402 | ! |
---|
| 403 | ! 0.1 input ! |
---|
| 404 | |
---|
| 405 | ! Domain size |
---|
| 406 | INTEGER(i_std), INTENT(in) :: npts_glo |
---|
| 407 | ! Size of cartesian grid |
---|
| 408 | INTEGER(i_std), INTENT(in) :: iim, jjm |
---|
| 409 | ! Longitudes on cartesian grid |
---|
[3447] | 410 | REAL(r_std), DIMENSION(iim,jjm), INTENT(in) :: grid_lon |
---|
[8] | 411 | ! Latitudes on cartesian grid |
---|
[3447] | 412 | REAL(r_std), DIMENSION(iim,jjm), INTENT(in) :: grid_lat |
---|
[8] | 413 | ! Index of land point on 2D map (in local position) |
---|
[3447] | 414 | INTEGER(i_std), DIMENSION(npts_glo), INTENT(in) :: kindex |
---|
| 415 | ! The fraction of continent in the grid box [0-1] |
---|
| 416 | REAL(r_std), DIMENSION(npts_glo), OPTIONAL, INTENT(in) :: contfrac_tmp |
---|
[8] | 417 | ! |
---|
| 418 | ! |
---|
| 419 | ! ========================================================================= |
---|
| 420 | |
---|
[2348] | 421 | IF ( printlev >= 4 ) WRITE(numout,*) 'Entering grid_stuff' |
---|
[8] | 422 | |
---|
| 423 | ! default resolution |
---|
[2348] | 424 | IF ( printlev >=2 ) WRITE(numout,*) 'grid stuff: default resolution (m): ',default_resolution |
---|
[8] | 425 | ! |
---|
| 426 | !- |
---|
| 427 | IF (is_root_prc) THEN |
---|
[3447] | 428 | ! |
---|
[5559] | 429 | CALL grid_topolylist(NbSegments, npts_glo, iim, jjm, grid_lon, grid_lat, kindex, & |
---|
[3447] | 430 | & global, corners_g, neighbours_g, headings_g, seglength_g, area_g, ilandindex, jlandindex) |
---|
| 431 | ! |
---|
| 432 | IF (PRESENT(contfrac_tmp)) THEN |
---|
| 433 | ! |
---|
| 434 | ! Transfer the contfrac into the array managed in this module. |
---|
| 435 | ! |
---|
| 436 | contfrac_g(:) = contfrac_tmp(:) |
---|
[8] | 437 | ENDIF |
---|
| 438 | ! |
---|
[3447] | 439 | ENDIF |
---|
| 440 | ! |
---|
| 441 | ! With this the description of the grid is complete and the information |
---|
| 442 | ! can be scattered to all processors. |
---|
| 443 | ! |
---|
| 444 | CALL grid_scatter() |
---|
| 445 | ! |
---|
[7259] | 446 | CALL bcast(neighbours_g) |
---|
| 447 | CALL bcast(resolution_g) |
---|
| 448 | ! |
---|
[3447] | 449 | IF ( printlev >= 3 ) WRITE(numout,*) 'Leaving grid_stuff' |
---|
| 450 | |
---|
| 451 | END SUBROUTINE grid_stuff |
---|
| 452 | !! |
---|
| 453 | !! ============================================================================================================================= |
---|
| 454 | !! SUBROUTINE: grid_topolylist |
---|
| 455 | !! |
---|
| 456 | !>\BRIEF This routine transforms a regular grid into a list of polygons which are defined by the following |
---|
| 457 | !! quantities : |
---|
| 458 | !! |
---|
| 459 | !! corners : the n vertices of the polugon in longitude and latitude |
---|
| 460 | !! neighbours : the neighbouring land grid box for each of the vertices and segments |
---|
| 461 | !! headings : the direction in which the neighbour is |
---|
| 462 | !! seglength : the lenght of each segment |
---|
| 463 | !! area : the area of the polygon |
---|
| 464 | !! ilindex, jlindex : provides the i,j coordinates of the mesh in the global grid. |
---|
| 465 | !! |
---|
| 466 | !! DESCRIPTION: |
---|
| 467 | !! |
---|
| 468 | !! \n |
---|
| 469 | !_ ============================================================================================================================== |
---|
| 470 | !! |
---|
[5559] | 471 | SUBROUTINE grid_topolylist(nbseg, nland, iim, jjm, grid_lon, grid_lat, kindex, & |
---|
[3447] | 472 | & globalg, corners_loc, neighbours_loc, headings_loc, seglength_loc, & |
---|
| 473 | & area_loc, ilindex_loc, jlindex_loc) |
---|
| 474 | ! |
---|
| 475 | ! 0 interface |
---|
| 476 | ! |
---|
| 477 | IMPLICIT NONE |
---|
| 478 | ! |
---|
| 479 | ! 0.1 input ! |
---|
| 480 | ! Number of segments for each polygon |
---|
| 481 | INTEGER(i_std), INTENT(in) :: nbseg |
---|
| 482 | ! Number of land points on the grid |
---|
| 483 | INTEGER(i_std), INTENT(in) :: nland |
---|
| 484 | ! Size of cartesian grid |
---|
| 485 | INTEGER(i_std), INTENT(in) :: iim, jjm |
---|
| 486 | ! Longitudes on cartesian grid |
---|
| 487 | REAL(r_std), DIMENSION(iim,jjm), INTENT(in) :: grid_lon |
---|
| 488 | ! Latitudes on cartesian grid |
---|
| 489 | REAL(r_std), DIMENSION(iim,jjm), INTENT(in) :: grid_lat |
---|
| 490 | ! Index of land point on 2D map (in local position) |
---|
| 491 | INTEGER(i_std), DIMENSION(nland), INTENT(in) :: kindex |
---|
| 492 | ! |
---|
| 493 | ! 0.2 Output |
---|
| 494 | ! |
---|
[4744] | 495 | LOGICAL, INTENT(inout) :: globalg |
---|
[3447] | 496 | ! |
---|
| 497 | REAL(r_std), DIMENSION(nland,nbseg,2), INTENT(out) :: corners_loc |
---|
| 498 | INTEGER(i_std), DIMENSION(nland,nbseg*2), INTENT(out) :: neighbours_loc |
---|
| 499 | REAL(r_std), DIMENSION(nland,nbseg*2), INTENT(out) :: headings_loc |
---|
| 500 | REAL(r_std), DIMENSION(nland,nbseg), INTENT(out) :: seglength_loc |
---|
| 501 | REAL(r_std), DIMENSION(nland), INTENT(out) :: area_loc |
---|
| 502 | INTEGER(i_std), DIMENSION(nland), INTENT(out) :: ilindex_loc, jlindex_loc |
---|
| 503 | ! |
---|
| 504 | ! 0.3 Local variables |
---|
| 505 | ! |
---|
| 506 | INTEGER(i_std) :: i, is, iss |
---|
| 507 | REAL(r_std), DIMENSION(nland,2) :: center |
---|
| 508 | REAL(r_std) :: maxdellon, mindellon, maxlon, minlon |
---|
| 509 | REAL(r_std), DIMENSION(nland,nbseg*2) :: lonpoly, latpoly |
---|
| 510 | ! |
---|
[5559] | 511 | IF ( grid_type == regular_lonlat ) THEN |
---|
[8] | 512 | ! |
---|
[3447] | 513 | ! If we are in regular Lon Lat, then we test just the longitude and see if we span 0-360deg. |
---|
| 514 | ! |
---|
| 515 | maxdellon=MAXVAL(ABS(grid_lon(1:iim-1,1)-grid_lon(2:iim,1))) |
---|
| 516 | mindellon=MINVAL(ABS(grid_lon(1:iim-1,1)-grid_lon(2:iim,1))) |
---|
| 517 | maxlon=MAXVAL(grid_lon(1:iim,1)) |
---|
| 518 | minlon=MINVAL(grid_lon(1:iim,1)) |
---|
| 519 | ! |
---|
| 520 | ! test if it could be a global grid on 0 -> 360 |
---|
| 521 | ! |
---|
| 522 | IF ( minlon > 0 .AND. maxlon > 180 ) THEN |
---|
[5419] | 523 | IF ( (minlon - (maxdellon + min_sechiba)) <= 0 .AND. (maxlon + (maxdellon + min_sechiba)) >= 360) THEN |
---|
[3447] | 524 | globalg = .TRUE. |
---|
[8] | 525 | ELSE |
---|
[3447] | 526 | globalg = .FALSE. |
---|
[8] | 527 | ENDIF |
---|
| 528 | ! |
---|
[3447] | 529 | ! Test if it could be a -180 to 180 grid |
---|
| 530 | ! |
---|
| 531 | ELSE IF ( minlon < 0 .AND. maxlon > 0 ) THEN |
---|
[5419] | 532 | IF ( (minlon - (maxdellon + min_sechiba)) <= -180 .AND. (maxlon + (maxdellon + min_sechiba)) >= 180) THEN |
---|
[3447] | 533 | globalg = .TRUE. |
---|
[8] | 534 | ELSE |
---|
[3447] | 535 | globalg = .FALSE. |
---|
[8] | 536 | ENDIF |
---|
[3447] | 537 | ! |
---|
| 538 | ! If neither condition is met then it cannot be global. |
---|
| 539 | ! |
---|
[8] | 540 | ELSE |
---|
[3447] | 541 | globalg = .FALSE. |
---|
[8] | 542 | ENDIF |
---|
[5559] | 543 | ELSE IF ( grid_type == regular_xy ) THEN |
---|
[8] | 544 | ! |
---|
[5559] | 545 | ! The hypothesis is that if we are using grid regular_xy then it is not a global grid |
---|
[8] | 546 | ! |
---|
[3447] | 547 | globalg = .FALSE. |
---|
| 548 | ELSE |
---|
| 549 | STOP "Unknown grid" |
---|
| 550 | ENDIF |
---|
| 551 | ! |
---|
| 552 | ! 2.0 Transform the grid into a list of polygones while keeping the neighbour relations |
---|
| 553 | ! between these polygones. |
---|
| 554 | ! |
---|
| 555 | ! Each polygone starts with a vertex and alternates vertices and mid-points of segments. |
---|
| 556 | ! |
---|
[3568] | 557 | IF (nland == 1) THEN |
---|
[4744] | 558 | CALL haversine_singlepointploy(iim, jjm, grid_lon, grid_lat, nland, kindex, globalg, & |
---|
[3568] | 559 | & nbseg, lonpoly, latpoly, center, & |
---|
| 560 | & neighbours_loc, ilindex_loc, jlindex_loc) |
---|
[5559] | 561 | ELSE IF ( grid_type == regular_lonlat ) THEN |
---|
[4744] | 562 | CALL haversine_reglatlontoploy(iim, jjm, grid_lon, grid_lat, nland, kindex, globalg, & |
---|
[3447] | 563 | & nbseg, lonpoly, latpoly, center, & |
---|
| 564 | & neighbours_loc, ilindex_loc, jlindex_loc) |
---|
[5559] | 565 | ELSE IF ( grid_type == regular_xy ) THEN |
---|
[3447] | 566 | CALL haversine_regxytoploy(iim, jjm, grid_lon, grid_lat, nland, kindex, proj_stack, & |
---|
| 567 | & nbseg, lonpoly, latpoly, center, & |
---|
| 568 | & neighbours_loc, ilindex_loc, jlindex_loc) |
---|
| 569 | ELSE |
---|
| 570 | STOP "Unknown grid" |
---|
| 571 | ENDIF |
---|
| 572 | ! |
---|
| 573 | ! Save the longitude and latitudes nbseg corners (=vertices) of the polygones |
---|
| 574 | ! |
---|
| 575 | DO i=1,nland |
---|
| 576 | DO is=1,nbseg |
---|
| 577 | iss=(is-1)*2+1 |
---|
| 578 | corners_loc(i,is,1) = lonpoly(i,iss) |
---|
| 579 | corners_loc(i,is,2) = latpoly(i,iss) |
---|
| 580 | ENDDO |
---|
| 581 | ENDDO |
---|
| 582 | ! |
---|
| 583 | ! Get the heading normal to the 4 segments and through the 4 corners. |
---|
| 584 | ! |
---|
| 585 | CALL haversine_polyheadings(nland, nbseg, lonpoly, latpoly, center, headings_loc) |
---|
| 586 | ! |
---|
| 587 | ! Order the points of the polygone in clockwise order Starting with the northern most |
---|
| 588 | ! |
---|
| 589 | CALL haversine_polysort(nland, nbseg, lonpoly, latpoly, headings_loc, neighbours_loc) |
---|
| 590 | ! |
---|
| 591 | ! Compute the segment length and area. |
---|
| 592 | ! For the RegLonLat we have specific calculations for seglength and area. |
---|
| 593 | ! For projected regular grids we use the great cicle assumption for the segments |
---|
| 594 | ! but the projected area. |
---|
| 595 | ! For unstructured grid we use the most general routines. |
---|
| 596 | ! |
---|
[5559] | 597 | IF ( grid_type == regular_lonlat ) THEN |
---|
[3447] | 598 | CALL haversine_laloseglen(nland, nbseg, lonpoly, latpoly, seglength_loc) |
---|
| 599 | CALL haversine_laloarea(nland, nbseg, seglength_loc, area_loc) |
---|
[5559] | 600 | ELSE IF ( grid_type == regular_xy ) THEN |
---|
[3447] | 601 | CALL haversine_polyseglen(nland, nbseg, lonpoly, latpoly, seglength_loc) |
---|
| 602 | CALL haversine_xyarea(nland, nbseg, ilindex_loc, jlindex_loc, dxwrf, dywrf, area_loc) |
---|
| 603 | ELSE |
---|
| 604 | CALL haversine_polyseglen(nland, nbseg, lonpoly, latpoly, seglength_loc) |
---|
| 605 | CALL haversine_polyarea(nland, nbseg, lonpoly, latpoly, area_loc) |
---|
| 606 | ENDIF |
---|
| 607 | |
---|
| 608 | END SUBROUTINE grid_topolylist |
---|
| 609 | !! |
---|
| 610 | !! |
---|
| 611 | !! |
---|
| 612 | !! ============================================================================================================================= |
---|
| 613 | !! SUBROUTINE: grid_scatter |
---|
| 614 | !! |
---|
| 615 | !>\BRIEF Scatter the grid information so that each processor knows the characteristics of the grid it works on. |
---|
| 616 | !! |
---|
| 617 | !! DESCRIPTION: |
---|
| 618 | !! |
---|
| 619 | !! |
---|
| 620 | !! The grid information has been computed for the entire grid on the root processor. Now we give each processor |
---|
| 621 | !! the information of the piece of the grid it works on. This concerns the following variables describing the grid : |
---|
| 622 | !! - area |
---|
| 623 | !! - resolution |
---|
| 624 | !! - neighbours |
---|
| 625 | !! - contfrac : fraction of continent |
---|
| 626 | !! |
---|
| 627 | !! Should ilandindex and jlandindex not b initialized, we catch-up here. This field is the same on all processors. |
---|
| 628 | !! |
---|
| 629 | !! TODO : |
---|
| 630 | !! This code should get the grid describing fields as arguments and then writem into the *_g variables on |
---|
| 631 | !! root_prc before scattering. This would allow to compute the grid characteristics in any subroutine |
---|
| 632 | !! fore calling grid_scatter. |
---|
| 633 | !! |
---|
| 634 | !! |
---|
| 635 | !! |
---|
| 636 | !! \n |
---|
| 637 | !_ ============================================================================================================================== |
---|
| 638 | !! |
---|
| 639 | !! |
---|
| 640 | SUBROUTINE grid_scatter() |
---|
| 641 | ! |
---|
| 642 | ! |
---|
| 643 | INTEGER(i_std) :: i, ip, jp |
---|
| 644 | ! |
---|
| 645 | IF ( MAXVAL(ilandindex) < 0 .AND. MAXVAL(jlandindex) < 0 ) THEN |
---|
| 646 | DO i = 1, nbp_glo |
---|
[8] | 647 | ! |
---|
| 648 | ! 1 find numbers of the latitude and longitude of each point |
---|
| 649 | ! |
---|
| 650 | |
---|
| 651 | ! index of latitude |
---|
[3447] | 652 | jp = INT( (index_g(i)-1) /iim_g ) + 1 |
---|
[8] | 653 | |
---|
| 654 | ! index of longitude |
---|
[3447] | 655 | ip = index_g(i) - ( jp-1 ) * iim_g |
---|
[8] | 656 | ! |
---|
[3447] | 657 | ! Save this information for usage in other modules. |
---|
[8] | 658 | ! |
---|
[3447] | 659 | ilandindex(i)=ip |
---|
| 660 | jlandindex(i)=jp |
---|
| 661 | ! |
---|
| 662 | ENDDO |
---|
| 663 | ENDIF |
---|
| 664 | ! |
---|
| 665 | CALL scatter(neighbours_g, neighbours) |
---|
| 666 | CALL scatter(contfrac_g, contfrac) |
---|
| 667 | CALL scatter(headings_g, headings) |
---|
| 668 | CALL scatter(seglength_g, seglength) |
---|
| 669 | CALL scatter(corners_g, corners) |
---|
| 670 | CALL scatter(area_g, area) |
---|
| 671 | ! |
---|
| 672 | ! TEMPORARY section for resolution |
---|
| 673 | ! |
---|
| 674 | IF ( is_root_prc) THEN |
---|
[5559] | 675 | IF ( grid_type == regular_lonlat .OR. grid_type == regular_xy ) THEN |
---|
[3447] | 676 | resolution_g(:,1) = (seglength_g(:,1)+seglength_g(:,3))/2.0 |
---|
| 677 | resolution_g(:,2) = (seglength_g(:,2)+seglength_g(:,4))/2.0 |
---|
[8] | 678 | ELSE |
---|
[3447] | 679 | CALL ipslerr(3, "grid_scatter", "unsupported grid type.",& |
---|
[5559] | 680 | & "grid_scatter can only be called for regular grids.",& |
---|
| 681 | & "") |
---|
[8] | 682 | ENDIF |
---|
[3447] | 683 | ENDIF |
---|
| 684 | CALL scatter(resolution_g, resolution) |
---|
[4744] | 685 | |
---|
| 686 | ! Copy variable global from root processor to all prossesors |
---|
| 687 | CALL bcast(global) |
---|
| 688 | |
---|
[3447] | 689 | ! |
---|
| 690 | ! |
---|
| 691 | IF ( printlev >=4 ) THEN |
---|
| 692 | WRITE(numout,*) 'grid_scatter > seglength = ', seglength(1,:) |
---|
| 693 | WRITE(numout,*) 'grid_scatter > neighbours = ', neighbours(1,:) |
---|
| 694 | WRITE(numout,*) 'grid_scatter > contfrac = ', contfrac(1) |
---|
| 695 | WRITE(numout,*) 'grid_scatter > area = ', area(1) |
---|
[4744] | 696 | WRITE(numout,*) 'grid_scatter > global = ', global |
---|
[3447] | 697 | ENDIF |
---|
| 698 | ! |
---|
| 699 | END SUBROUTINE grid_scatter |
---|
| 700 | !! |
---|
| 701 | !! |
---|
| 702 | !! ============================================================================================================================= |
---|
| 703 | !! SUBROUTINE: grid_initproj |
---|
| 704 | !! |
---|
| 705 | !>\BRIEF Routine to initialise the projection |
---|
| 706 | !! |
---|
| 707 | !! DESCRIPTION: |
---|
| 708 | !! |
---|
| 709 | !! |
---|
| 710 | !! This subroutine is called by the routine whichs ets-up th grid on which ORCHIDEE is to run. |
---|
| 711 | !! The aim is to set-upu the projection so that all the grid variables needed by ORCHIDEE can |
---|
| 712 | !! be computed in grid_stuff_regxy |
---|
| 713 | !! |
---|
| 714 | !! \n |
---|
| 715 | !_ ============================================================================================================================== |
---|
| 716 | !! |
---|
| 717 | !! |
---|
| 718 | SUBROUTINE grid_initproj (fid, iim, jjm) |
---|
| 719 | ! |
---|
| 720 | ! |
---|
| 721 | ! 0 interface |
---|
| 722 | ! |
---|
| 723 | IMPLICIT NONE |
---|
| 724 | ! |
---|
| 725 | ! 0.1 input ! |
---|
| 726 | ! |
---|
| 727 | ! Domain size |
---|
| 728 | INTEGER(i_std), INTENT(in) :: fid |
---|
| 729 | INTEGER(i_std), INTENT(in) :: iim, jjm |
---|
| 730 | ! |
---|
| 731 | ! 0.2 Local variables |
---|
| 732 | ! |
---|
| 733 | INTEGER(i_std) :: current_proj, idom, iret, lonid, latid, numLons, numLats |
---|
| 734 | INTEGER, DIMENSION(nf90_max_var_dims) :: dimIDs |
---|
| 735 | REAL(r_std) :: user_stand_lon, user_truelat1, user_truelat2, user_dxkm, user_dykm |
---|
| 736 | REAL(r_std) :: user_dlat, user_dlon, user_known_x, user_known_y, user_known_lat, user_known_lon |
---|
| 737 | REAL(r_std), DIMENSION(16) :: corner_lons, corner_lats |
---|
| 738 | ! |
---|
| 739 | INTEGER(i_std) :: iv, i, j |
---|
| 740 | CHARACTER(LEN=20) :: varname |
---|
| 741 | REAL(r_std) :: dx, dy, dtx, dty, coslat |
---|
| 742 | REAL(r_std), ALLOCATABLE, DIMENSION (:) :: LON, LAT |
---|
| 743 | REAL(r_std), ALLOCATABLE, DIMENSION (:,:) :: mapfac_x, mapfac_y |
---|
| 744 | ! |
---|
| 745 | ! |
---|
| 746 | ! Only one domain is possible for the moment |
---|
| 747 | ! |
---|
| 748 | idom=1 |
---|
| 749 | CALL map_init(proj_stack(idom)) |
---|
| 750 | ! |
---|
| 751 | ! Does ORCHIDEE have the same Earth Radius as the map projection ? |
---|
| 752 | ! |
---|
| 753 | IF ( ABS(R_Earth-EARTH_RADIUS_M) > 0.1 ) THEN |
---|
| 754 | WRITE(*,*) "Earth Radius in WRF : ", EARTH_RADIUS_M |
---|
| 755 | WRITE(*,*) "Earth Radius in ORCHIDEE : ", R_Earth |
---|
| 756 | CALL ipslerr (3,'grid_initproj','The Earth radius is not the same in the projection module and ORCHIDEE',& |
---|
| 757 | & " ", " ") |
---|
| 758 | ENDIF |
---|
| 759 | ! |
---|
| 760 | ! Get parameters of the projection from the netCDF file |
---|
| 761 | ! |
---|
| 762 | iret = NF90_GET_ATT(fid, NF90_GLOBAL, "MAP_PROJ", current_proj) |
---|
| 763 | ! |
---|
| 764 | iret = NF90_GET_ATT(fid, NF90_GLOBAL, "STAND_LON", user_stand_lon) |
---|
| 765 | iret = NF90_GET_ATT(fid, NF90_GLOBAL, "TRUELAT1", user_truelat1) |
---|
| 766 | iret = NF90_GET_ATT(fid, NF90_GLOBAL, "TRUELAT2", user_truelat2) |
---|
| 767 | ! |
---|
| 768 | iret = NF90_GET_ATT(fid, NF90_GLOBAL, "DX", user_dxkm) |
---|
| 769 | iret = NF90_GET_ATT(fid, NF90_GLOBAL, "DY", user_dykm) |
---|
| 770 | user_dlat = undef |
---|
| 771 | user_dlon = undef |
---|
| 772 | ! |
---|
| 773 | IF ( current_proj == PROJ_LATLON ) THEN |
---|
| 774 | ! |
---|
| 775 | iret = NF90_inq_VARID(fid, "XLONG_M",lonid) |
---|
| 776 | iret = NF90_INQUIRE_VARIABLE(fid, lonid, dimids = dimIDs) |
---|
| 777 | iret = NF90_INQUIRE_DIMENSION(fid, dimIDs(1), len = numLons) |
---|
| 778 | iret = NF90_INQUIRE_DIMENSION(fid, dimIDs(2), len = numLats) |
---|
| 779 | ALLOCATE(LON(numLons)) |
---|
| 780 | iret = NF90_GET_VAR(fid, lonid, LON(:), start = (/ 1, 1, 1 /), count = (/ numLons, 1, 1 /)) |
---|
| 781 | |
---|
| 782 | iret = NF90_inq_VARID(fid, "XLAT_M",latid) |
---|
| 783 | ALLOCATE(LAT(numLats)) |
---|
| 784 | iret = NF90_GET_VAR(fid, latid, LAT(:), start = (/ 1, 1, 1 /), count = (/ 1, numLats, 1 /)) |
---|
| 785 | |
---|
| 786 | user_dlon = (LON(numLons) - LON(1)) / (numLons - 1) |
---|
| 787 | user_dlat = (LAT(numLats) - LAT(1)) / (numLats - 1) |
---|
| 788 | |
---|
| 789 | DEALLOCATE(LON,LAT) |
---|
[8] | 790 | |
---|
[3447] | 791 | ENDIF |
---|
| 792 | ! Unable to know from where to get the information |
---|
| 793 | user_known_x = 1 |
---|
| 794 | user_known_y = 1 |
---|
| 795 | ! |
---|
| 796 | iret = NF90_GET_ATT(fid, NF90_GLOBAL, "corner_lats", corner_lats) |
---|
| 797 | iret = NF90_GET_ATT(fid, NF90_GLOBAL, "corner_lons", corner_lons) |
---|
| 798 | user_known_lat = corner_lats(1) |
---|
| 799 | user_known_lon = corner_lons(1) |
---|
| 800 | ! |
---|
| 801 | ! Read mapfactor, land mask and orography |
---|
| 802 | ! |
---|
| 803 | ! |
---|
| 804 | ! Allocation |
---|
| 805 | ! |
---|
| 806 | ALLOCATE(mapfac_x(iim,jjm)) |
---|
| 807 | ALLOCATE(mapfac_y(iim,jjm)) |
---|
| 808 | ALLOCATE(dxwrf(iim,jjm)) |
---|
| 809 | ALLOCATE(dywrf(iim,jjm)) |
---|
| 810 | ! |
---|
| 811 | varname = "MAPFAC_MX" |
---|
| 812 | iret = NF90_INQ_VARID (fid, varname, iv) |
---|
| 813 | IF (iret /= NF90_NOERR) THEN |
---|
| 814 | CALL ipslerr (3,'WRFdomain_Read',"Could not find variable ", varname," ") |
---|
| 815 | ELSE |
---|
| 816 | iret = NF90_GET_VAR (fid,iv,mapfac_x) |
---|
| 817 | ENDIF |
---|
| 818 | varname = "MAPFAC_MY" |
---|
| 819 | iret = NF90_INQ_VARID (fid, varname, iv) |
---|
| 820 | IF (iret /= NF90_NOERR) THEN |
---|
| 821 | CALL ipslerr (3,'WRFdomain_Read',"Could not find variable ", varname," ") |
---|
| 822 | ELSE |
---|
| 823 | iret = NF90_GET_VAR (fid,iv,mapfac_y) |
---|
| 824 | ENDIF |
---|
| 825 | ! |
---|
| 826 | ! Initilize the projection |
---|
| 827 | ! |
---|
| 828 | if (current_proj == PROJ_LATLON) then |
---|
| 829 | call map_set(current_proj, proj_stack(idom), & |
---|
| 830 | lat1=user_known_lat, & |
---|
| 831 | lon1=user_known_lon, & |
---|
| 832 | knowni=user_known_x, & |
---|
| 833 | knownj=user_known_y, & |
---|
| 834 | latinc=user_dlat, & |
---|
| 835 | loninc=user_dlon, & |
---|
| 836 | r_earth=R_Earth) |
---|
| 837 | |
---|
| 838 | else if (current_proj == PROJ_MERC) then |
---|
| 839 | call map_set(current_proj, proj_stack(idom), & |
---|
| 840 | truelat1=user_truelat1, & |
---|
| 841 | lat1=user_known_lat, & |
---|
| 842 | lon1=user_known_lon, & |
---|
| 843 | knowni=user_known_x, & |
---|
| 844 | knownj=user_known_y, & |
---|
| 845 | dx=user_dxkm, & |
---|
| 846 | r_earth=R_Earth) |
---|
| 847 | |
---|
| 848 | else if (current_proj == PROJ_CYL) then |
---|
| 849 | call ipslerr(3,"grid_initproj",'Should not have PROJ_CYL as projection for',& |
---|
| 850 | 'source data in push_source_projection()', " ") |
---|
| 851 | |
---|
| 852 | else if (current_proj == PROJ_CASSINI) then |
---|
| 853 | call ipslerr(3,"grid_initproj",'Should not have PROJ_CASSINI as projection for', & |
---|
| 854 | 'source data in push_source_projection()', " ") |
---|
| 855 | |
---|
| 856 | else if (current_proj == PROJ_LC) then |
---|
| 857 | call map_set(current_proj, proj_stack(idom), & |
---|
| 858 | truelat1=user_truelat1, & |
---|
| 859 | truelat2=user_truelat2, & |
---|
| 860 | stdlon=user_stand_lon, & |
---|
| 861 | lat1=user_known_lat, & |
---|
| 862 | lon1=user_known_lon, & |
---|
| 863 | knowni=user_known_x, & |
---|
| 864 | knownj=user_known_y, & |
---|
| 865 | dx=user_dxkm, & |
---|
| 866 | r_earth=R_Earth) |
---|
| 867 | |
---|
| 868 | else if (current_proj == PROJ_ALBERS_NAD83) then |
---|
| 869 | call map_set(current_proj, proj_stack(idom), & |
---|
| 870 | truelat1=user_truelat1, & |
---|
| 871 | truelat2=user_truelat2, & |
---|
| 872 | stdlon=user_stand_lon, & |
---|
| 873 | lat1=user_known_lat, & |
---|
| 874 | lon1=user_known_lon, & |
---|
| 875 | knowni=user_known_x, & |
---|
| 876 | knownj=user_known_y, & |
---|
| 877 | dx=user_dxkm, & |
---|
| 878 | r_earth=R_Earth) |
---|
| 879 | |
---|
| 880 | else if (current_proj == PROJ_PS) then |
---|
| 881 | call map_set(current_proj, proj_stack(idom), & |
---|
| 882 | truelat1=user_truelat1, & |
---|
| 883 | stdlon=user_stand_lon, & |
---|
| 884 | lat1=user_known_lat, & |
---|
| 885 | lon1=user_known_lon, & |
---|
| 886 | knowni=user_known_x, & |
---|
| 887 | knownj=user_known_y, & |
---|
| 888 | dx=user_dxkm, & |
---|
| 889 | r_earth=R_Earth) |
---|
| 890 | |
---|
| 891 | else if (current_proj == PROJ_PS_WGS84) then |
---|
| 892 | call map_set(current_proj, proj_stack(idom), & |
---|
| 893 | truelat1=user_truelat1, & |
---|
| 894 | stdlon=user_stand_lon, & |
---|
| 895 | lat1=user_known_lat, & |
---|
| 896 | lon1=user_known_lon, & |
---|
| 897 | knowni=user_known_x, & |
---|
| 898 | knownj=user_known_y, & |
---|
| 899 | dx=user_dxkm, & |
---|
| 900 | r_earth=R_Earth) |
---|
| 901 | |
---|
| 902 | else if (current_proj == PROJ_GAUSS) then |
---|
| 903 | call map_set(current_proj, proj_stack(idom), & |
---|
| 904 | lat1=user_known_lat, & |
---|
| 905 | lon1=user_known_lon, & |
---|
| 906 | nlat=nint(user_dlat), & |
---|
| 907 | loninc=user_dlon, & |
---|
| 908 | r_earth=R_Earth) |
---|
| 909 | |
---|
| 910 | else if (current_proj == PROJ_ROTLL) then |
---|
| 911 | call ipslerr(3 ,"grid_initproj",'Should not have PROJ_ROTLL as projection for', & |
---|
| 912 | 'source data in push_source_projection() as not yet implemented', '') |
---|
| 913 | end if |
---|
| 914 | ! |
---|
| 915 | ! Transform the mapfactors into dx and dy to be used for the description of the polygons and |
---|
| 916 | ! interpolations. |
---|
| 917 | ! |
---|
| 918 | DO i=1,iim |
---|
| 919 | DO j=1,jjm |
---|
[8] | 920 | ! |
---|
[3447] | 921 | IF (proj_stack(idom)%code /= PROJ_LATLON ) THEN |
---|
| 922 | dx = proj_stack(idom)%dx |
---|
| 923 | ! Some projections in WRF do not store dy, in that case dy=dx. |
---|
| 924 | IF ( proj_stack(idom)%dy > 0 ) THEN |
---|
| 925 | dy = proj_stack(idom)%dy |
---|
[8] | 926 | ELSE |
---|
[3447] | 927 | dy = proj_stack(idom)%dx |
---|
[8] | 928 | ENDIF |
---|
[3447] | 929 | dxwrf(i,j) = dx/mapfac_x(i,j) |
---|
| 930 | dywrf(i,j) = dy/mapfac_y(i,j) |
---|
[8] | 931 | ELSE |
---|
[3447] | 932 | ! |
---|
| 933 | ! The LatLon projection is also a special case as here it is not the dx and dy |
---|
| 934 | ! which are stored in the projection file but the increments in Lon and Lat. |
---|
| 935 | ! |
---|
| 936 | dtx = proj_stack(idom)%loninc |
---|
| 937 | dty = proj_stack(idom)%latinc |
---|
| 938 | coslat = COS(lat(j) * pi/180. ) |
---|
| 939 | dxwrf(i,j) = dtx * pi/180. * R_Earth * coslat |
---|
| 940 | dywrf(i,j) = dty * pi/180. * R_Earth |
---|
| 941 | ! |
---|
[8] | 942 | ENDIF |
---|
| 943 | ! |
---|
[3447] | 944 | ENDDO |
---|
| 945 | ENDDO |
---|
| 946 | ! |
---|
| 947 | END SUBROUTINE grid_initproj |
---|
| 948 | ! |
---|
| 949 | ! |
---|
| 950 | ! |
---|
| 951 | !========================================================================================= |
---|
| 952 | ! |
---|
| 953 | SUBROUTINE grid_tolola_scal (ri, rj, lon, lat) |
---|
| 954 | ! |
---|
| 955 | ! |
---|
| 956 | ! Argument |
---|
| 957 | REAL(r_std), INTENT(in) :: ri, rj |
---|
| 958 | REAL(r_std), INTENT(out) :: lon, lat |
---|
| 959 | ! |
---|
| 960 | ! |
---|
| 961 | IF ( proj_stack(1)%code < undef_int ) THEN |
---|
| 962 | ! |
---|
| 963 | CALL ij_to_latlon(proj_stack(1), ri, rj, lat, lon) |
---|
| 964 | ! |
---|
| 965 | ELSE |
---|
| 966 | CALL ipslerr(3, "grid_tolola_scal", "Projection not initilized"," "," ") |
---|
| 967 | ENDIF |
---|
| 968 | ! |
---|
| 969 | END SUBROUTINE grid_tolola_scal |
---|
| 970 | ! |
---|
| 971 | !========================================================================================= |
---|
| 972 | ! |
---|
| 973 | SUBROUTINE grid_tolola_1d (ri, rj, lon, lat) |
---|
| 974 | ! |
---|
| 975 | ! |
---|
| 976 | ! Argument |
---|
| 977 | REAL(r_std), INTENT(in), DIMENSION(:) :: ri, rj |
---|
| 978 | REAL(r_std), INTENT(out), DIMENSION(:) :: lon, lat |
---|
| 979 | ! |
---|
| 980 | ! Local |
---|
| 981 | INTEGER :: i, imax |
---|
| 982 | ! |
---|
| 983 | imax=SIZE(lon) |
---|
| 984 | ! |
---|
| 985 | IF ( proj_stack(1)%code < undef_int ) THEN |
---|
| 986 | DO i=1,imax |
---|
[8] | 987 | ! |
---|
[3447] | 988 | CALL ij_to_latlon(proj_stack(1), ri(i), rj(i), lat(i), lon(i)) |
---|
[8] | 989 | ! |
---|
[3447] | 990 | ENDDO |
---|
| 991 | ELSE |
---|
| 992 | CALL ipslerr(3, "grid_tolola_1d", "Projection not initilized"," "," ") |
---|
| 993 | ENDIF |
---|
| 994 | ! |
---|
| 995 | END SUBROUTINE grid_tolola_1d |
---|
| 996 | ! |
---|
| 997 | !========================================================================================= |
---|
| 998 | ! |
---|
| 999 | SUBROUTINE grid_tolola_2d (ri, rj, lon, lat) |
---|
| 1000 | ! |
---|
| 1001 | ! |
---|
| 1002 | ! Argument |
---|
| 1003 | REAL(r_std), INTENT(in), DIMENSION(:,:) :: ri, rj |
---|
| 1004 | REAL(r_std), INTENT(out), DIMENSION(:,:) :: lon, lat |
---|
| 1005 | ! |
---|
| 1006 | ! Local |
---|
| 1007 | INTEGER :: i, imax, j, jmax |
---|
| 1008 | ! |
---|
| 1009 | imax=SIZE(lon,DIM=1) |
---|
| 1010 | jmax=SIZE(lon,DIM=2) |
---|
| 1011 | ! |
---|
| 1012 | IF ( proj_stack(1)%code < undef_int ) THEN |
---|
| 1013 | DO i=1,imax |
---|
| 1014 | DO j=1,jmax |
---|
| 1015 | ! |
---|
| 1016 | CALL ij_to_latlon(proj_stack(1), ri(i,j), rj(i,j), lat(i,j), lon(i,j)) |
---|
| 1017 | ! |
---|
| 1018 | ENDDO |
---|
| 1019 | ENDDO |
---|
| 1020 | ELSE |
---|
| 1021 | CALL ipslerr(3, "grid_tolola_2d", "Projection not initilized"," "," ") |
---|
| 1022 | ENDIF |
---|
| 1023 | ! |
---|
| 1024 | END SUBROUTINE grid_tolola_2d |
---|
| 1025 | ! |
---|
| 1026 | !========================================================================================= |
---|
| 1027 | ! |
---|
| 1028 | SUBROUTINE grid_toij_scal (lon, lat, ri, rj) |
---|
| 1029 | ! |
---|
| 1030 | ! |
---|
| 1031 | ! Argument |
---|
| 1032 | REAL(r_std), INTENT(in) :: lon, lat |
---|
| 1033 | REAL(r_std), INTENT(out) :: ri, rj |
---|
| 1034 | ! |
---|
| 1035 | ! |
---|
| 1036 | IF ( proj_stack(1)%code < undef_int ) THEN |
---|
| 1037 | ! |
---|
| 1038 | CALL latlon_to_ij(proj_stack(1), lat, lon, ri, rj) |
---|
| 1039 | ! |
---|
| 1040 | ELSE |
---|
| 1041 | CALL ipslerr(3, "grid_toij_scal", "Projection not initilized"," "," ") |
---|
| 1042 | ENDIF |
---|
| 1043 | ! |
---|
| 1044 | END SUBROUTINE grid_toij_scal |
---|
| 1045 | ! |
---|
| 1046 | !========================================================================================= |
---|
| 1047 | ! |
---|
| 1048 | SUBROUTINE grid_toij_1d (lon, lat, ri, rj) |
---|
| 1049 | ! |
---|
| 1050 | ! |
---|
| 1051 | ! Argument |
---|
| 1052 | REAL(r_std), INTENT(in), DIMENSION(:) :: lon, lat |
---|
| 1053 | REAL(r_std), INTENT(out), DIMENSION(:) :: ri, rj |
---|
| 1054 | ! |
---|
| 1055 | ! Local |
---|
| 1056 | INTEGER :: i, imax |
---|
| 1057 | ! |
---|
| 1058 | imax=SIZE(lon) |
---|
| 1059 | ! |
---|
| 1060 | IF ( proj_stack(1)%code < undef_int ) THEN |
---|
| 1061 | DO i=1,imax |
---|
[8] | 1062 | ! |
---|
[3447] | 1063 | CALL latlon_to_ij(proj_stack(1), lat(i), lon(i), ri(i), rj(i)) |
---|
[8] | 1064 | ! |
---|
| 1065 | ENDDO |
---|
[3447] | 1066 | ELSE |
---|
| 1067 | CALL ipslerr(3, "grid_toij_1d", "Projection not initilized"," "," ") |
---|
[8] | 1068 | ENDIF |
---|
[3447] | 1069 | ! |
---|
| 1070 | END SUBROUTINE grid_toij_1d |
---|
| 1071 | ! |
---|
| 1072 | !========================================================================================= |
---|
| 1073 | ! |
---|
| 1074 | SUBROUTINE grid_toij_2d (lon, lat, ri, rj) |
---|
| 1075 | ! |
---|
| 1076 | ! |
---|
| 1077 | ! Argument |
---|
| 1078 | REAL(r_std), INTENT(in), DIMENSION(:,:) :: lon, lat |
---|
| 1079 | REAL(r_std), INTENT(out), DIMENSION(:,:) :: ri, rj |
---|
| 1080 | ! |
---|
| 1081 | ! Local |
---|
| 1082 | INTEGER :: i, imax, j, jmax |
---|
| 1083 | ! |
---|
| 1084 | imax=SIZE(lon,DIM=1) |
---|
| 1085 | jmax=SIZE(lon,DIM=2) |
---|
| 1086 | ! |
---|
| 1087 | IF ( proj_stack(1)%code < undef_int ) THEN |
---|
| 1088 | DO i=1,imax |
---|
| 1089 | DO j=1,jmax |
---|
| 1090 | ! |
---|
| 1091 | CALL latlon_to_ij(proj_stack(1), lat(i,j), lon(i,j), ri(i,j), rj(i,j)) |
---|
| 1092 | ! |
---|
| 1093 | ENDDO |
---|
| 1094 | ENDDO |
---|
| 1095 | ELSE |
---|
| 1096 | CALL ipslerr(3, "grid_toij_2d", "Projection not initilized"," "," ") |
---|
| 1097 | ENDIF |
---|
| 1098 | ! |
---|
| 1099 | END SUBROUTINE grid_toij_2d |
---|
| 1100 | ! |
---|
| 1101 | ! |
---|
| 1102 | !========================================================================================= |
---|
| 1103 | ! |
---|
| 1104 | ! |
---|
[8] | 1105 | END MODULE grid |
---|