[2136] | 1 | !************************************************************************ |
---|
| 2 | ! Fortran 95 OPA Nesting tools * |
---|
| 3 | ! * |
---|
| 4 | ! Copyright (C) 2005 Florian Lemarié (Florian.Lemarie@imag.fr) * |
---|
| 5 | ! Laurent Debreu (Laurent.Debreu@imag.fr) * |
---|
| 6 | !************************************************************************ |
---|
| 7 | ! |
---|
| 8 | PROGRAM create_bathy |
---|
| 9 | ! |
---|
| 10 | USE NETCDF |
---|
| 11 | USE bilinear_interp |
---|
| 12 | USE agrif_readwrite |
---|
| 13 | USE agrif_partial_steps |
---|
| 14 | USE agrif_connect_topo |
---|
| 15 | ! |
---|
| 16 | IMPLICIT NONE |
---|
| 17 | ! |
---|
| 18 | !************************************************************************ |
---|
| 19 | ! * |
---|
[2455] | 20 | ! PROGRAM CREATE_BATHY * |
---|
[2136] | 21 | ! * |
---|
| 22 | ! program to implement bathymetry interpolation to generate * |
---|
| 23 | ! child grid bathymetry file * |
---|
| 24 | ! * |
---|
| 25 | ! various options : * |
---|
| 26 | ! * |
---|
| 27 | ! 1- Interpolation directly from parent bathymetry file (z-coord) * |
---|
[2455] | 28 | ! 2- Use new topo file in meters (for example etopo) * |
---|
[2136] | 29 | ! * |
---|
| 30 | ! vertical coordinates permitted : z-coord and partial steps * |
---|
| 31 | ! sigma coordinates is not yet implemented * |
---|
| 32 | ! * |
---|
| 33 | !Interpolation is carried out using bilinear interpolation * |
---|
| 34 | !routine from SCRIP package or median average * |
---|
| 35 | ! * |
---|
| 36 | !http://climate.lanl.gov/Software/SCRIP/ * |
---|
| 37 | !************************************************************************ |
---|
| 38 | ! |
---|
| 39 | ! variables declaration |
---|
| 40 | ! |
---|
| 41 | CHARACTER(len=80) :: namelistname |
---|
| 42 | CHARACTER*100 :: Childmeter_file,Childlevel_file,Child_coordinates,child_ps |
---|
| 43 | LOGICAL,DIMENSION(:,:),POINTER :: masksrc => NULL() |
---|
| 44 | LOGICAL :: identical_grids |
---|
| 45 | INTEGER,DIMENSION(:,:),ALLOCATABLE ::mask_oce,trouble_points |
---|
| 46 | INTEGER :: i,j,num_links,nb,nbadd,status,narg,iargc |
---|
| 47 | INTEGER,DIMENSION(:),POINTER :: src_add,dst_add => NULL() |
---|
| 48 | INTEGER :: numlatfine,numlonfine,numlat,numlon,pos,pos2 |
---|
| 49 | REAL*8,DIMENSION(:,:),POINTER :: matrix,interpdata => NULL() |
---|
| 50 | REAL*8, DIMENSION(:,:),POINTER :: bathy_fin_constant => NULL() |
---|
| 51 | REAL*8,DIMENSION(:,:),ALLOCATABLE :: bathy_test,vardep,glamhr,gphihr |
---|
| 52 | REAL*8,DIMENSION(:),ALLOCATABLE :: vardep1d |
---|
| 53 | REAL*8, DIMENSION(:,:),POINTER :: gdepw_ps_interp => NULL() |
---|
| 54 | REAL*8, DIMENSION(:,:),POINTER :: save_gdepw,rx,ry,maskedtopo |
---|
[9149] | 55 | REAL*8, DIMENSION(:,:),POINTER :: wtab => NULL() |
---|
[2136] | 56 | REAL*8 :: Cell_lonmin,Cell_lonmax,Cell_latmin,Cell_latmax,wghts |
---|
[8641] | 57 | LOGICAL :: Pacifique=.FALSE. |
---|
[2136] | 58 | INTEGER :: boundary,xpos,ypos,iimin,iimax,jjmax,jjmin |
---|
| 59 | INTEGER :: nbloops,nxhr,nyhr,ji,jj,nbiter,nbloopmax |
---|
| 60 | INTEGER :: ipt,jpt,iloc,jloc |
---|
| 61 | INTEGER, DIMENSION(2) :: i_min,i_max,j_min,j_max |
---|
| 62 | |
---|
| 63 | TYPE(Coordinates) :: G0,G1 |
---|
| 64 | ! |
---|
| 65 | narg = iargc() |
---|
| 66 | IF (narg == 0) THEN |
---|
| 67 | namelistname = 'namelist.input' |
---|
| 68 | ELSE |
---|
| 69 | CALL getarg(1,namelistname) |
---|
| 70 | ENDIF |
---|
| 71 | ! |
---|
| 72 | ! read input file (namelist.input) |
---|
| 73 | ! |
---|
| 74 | CALL read_namelist(namelistname) |
---|
| 75 | ! |
---|
| 76 | ! define names of child grid files |
---|
| 77 | ! |
---|
| 78 | CALL set_child_name(parent_coordinate_file,Child_coordinates) |
---|
[2455] | 79 | IF( TRIM(parent_meshmask_file) .NE. '/NULL' ) & |
---|
| 80 | CALL set_child_name(parent_meshmask_file,Childlevel_file) |
---|
[2136] | 81 | ! |
---|
| 82 | ! |
---|
| 83 | ! |
---|
| 84 | ! |
---|
| 85 | ! |
---|
| 86 | !------------------------------------------------------------------ |
---|
| 87 | ! ****First option : no new topo file / no partial steps |
---|
| 88 | ! interpolate levels directly from parent file |
---|
| 89 | !------------------------------------------------------------------ |
---|
| 90 | ! |
---|
| 91 | ! |
---|
| 92 | ! |
---|
| 93 | ! |
---|
| 94 | ! |
---|
| 95 | ! |
---|
| 96 | IF(.NOT.new_topo .AND. .NOT.partial_steps) THEN |
---|
| 97 | ! |
---|
| 98 | WRITE(*,*) 'First option' |
---|
| 99 | !read coarse grid bathymetry and coordinates file |
---|
| 100 | ! |
---|
| 101 | WRITE(*,*) 'No new topo file ...' |
---|
| 102 | status = Read_Coordinates(TRIM(parent_coordinate_file),G0) |
---|
[2455] | 103 | status = Read_bathy_level(TRIM(parent_meshmask_file),G0) |
---|
[2136] | 104 | ! |
---|
| 105 | IF( imax > SIZE(G0%glamt,1) .OR. jmax > SIZE(G0%glamt,2) .OR. & |
---|
| 106 | imax <= imin .OR. jmax <= jmin ) THEN |
---|
| 107 | WRITE(*,*) 'ERROR ***** bad child grid definition ...' |
---|
| 108 | WRITE(*,*) 'please check imin,jmin,imax,jmax,jpizoom,jpjzoom values' |
---|
| 109 | WRITE(*,*) ' ' |
---|
| 110 | STOP |
---|
| 111 | ENDIF |
---|
| 112 | ! |
---|
| 113 | !read fine grid coordinates file |
---|
| 114 | ! |
---|
| 115 | status = Read_Coordinates(TRIM(Child_coordinates),G1,pacifique) |
---|
| 116 | ! |
---|
| 117 | IF( SIZE(G1%nav_lon,1) .NE. nxfin .OR. SIZE(G1%nav_lon,2) .NE. nyfin ) THEN |
---|
| 118 | ! |
---|
| 119 | WRITE(*,*) 'ERROR ***** bad child coordinates file ...' |
---|
| 120 | WRITE(*,*) ' ' |
---|
| 121 | WRITE(*,*) 'please check that child coordinates file ' |
---|
| 122 | WRITE(*,*) 'has been created with the same namelist ' |
---|
| 123 | WRITE(*,*) ' ' |
---|
| 124 | STOP |
---|
| 125 | ! |
---|
| 126 | ENDIF |
---|
| 127 | ! |
---|
| 128 | numlat = SIZE(G0%nav_lat,2) |
---|
| 129 | numlon = SIZE(G0%nav_lat,1) |
---|
| 130 | numlatfine = SIZE(G1%nav_lat,2) |
---|
| 131 | numlonfine = SIZE(G1%nav_lat,1) |
---|
| 132 | ! |
---|
| 133 | ALLOCATE(masksrc(numlon,numlat)) |
---|
| 134 | ! |
---|
| 135 | ! create logical array masksrc |
---|
| 136 | ! |
---|
| 137 | WHERE(G0%bathy_level.LE.0) |
---|
| 138 | masksrc = .FALSE. |
---|
| 139 | ELSEWHERE |
---|
| 140 | masksrc = .TRUE. |
---|
| 141 | END WHERE |
---|
| 142 | |
---|
| 143 | IF ( Pacifique ) THEN |
---|
| 144 | WHERE(G0%nav_lon < 0.001) |
---|
| 145 | G0%nav_lon = G0%nav_lon + 360. |
---|
| 146 | END WHERE |
---|
| 147 | ENDIF |
---|
| 148 | |
---|
| 149 | |
---|
| 150 | !----------------- |
---|
| 151 | ! compute remapping matrix thanks to SCRIP package |
---|
| 152 | ! |
---|
| 153 | ! remapping process |
---|
| 154 | ! |
---|
| 155 | ALLOCATE(G1%bathy_meter(nxfin,nyfin)) |
---|
| 156 | CALL levels_to_meter(G0) |
---|
| 157 | ! |
---|
| 158 | ! Call levels_to_meter(G1) |
---|
| 159 | ! |
---|
| 160 | CALL get_remap_matrix(G0%nav_lat,G1%nav_lat, & |
---|
| 161 | G0%nav_lon,G1%nav_lon, & |
---|
| 162 | masksrc,matrix,src_add,dst_add) |
---|
| 163 | CALL make_remap(G0%bathy_meter,G1%bathy_meter,nxfin,nyfin, & |
---|
| 164 | matrix,src_add,dst_add) |
---|
| 165 | ! |
---|
| 166 | ! |
---|
| 167 | DEALLOCATE(masksrc) |
---|
| 168 | !----------------- |
---|
| 169 | ! |
---|
| 170 | ! |
---|
| 171 | ! compute constant bathymetry for Parent-Child bathymetry connection |
---|
| 172 | ! |
---|
| 173 | CALL init_constant_bathy(G0%bathy_meter,bathy_fin_constant) |
---|
| 174 | ! |
---|
[9166] | 175 | boundary = connectionsize*irafx + nbghostcellsfine + 1 |
---|
[2136] | 176 | ! |
---|
| 177 | ! connection carried out by copying parent grid values for the fine points |
---|
| 178 | ! corresponding to 3 coarse grid cells at the boundaries |
---|
| 179 | ! |
---|
| 180 | G1%bathy_meter(1:boundary,:) = bathy_fin_constant(1:boundary,:) |
---|
| 181 | G1%bathy_meter(:,1:boundary) = bathy_fin_constant(:,1:boundary) |
---|
| 182 | G1%bathy_meter(nxfin-boundary+1:nxfin,:) = bathy_fin_constant(nxfin-boundary+1:nxfin,:) |
---|
| 183 | G1%bathy_meter(:,nyfin-boundary+1:nyfin) = bathy_fin_constant(:,nyfin-boundary+1:nyfin) |
---|
| 184 | ! |
---|
| 185 | CALL smooth_topo(G1%bathy_meter(boundary:nxfin-boundary+1,boundary:nyfin-boundary+1),nbiter) |
---|
| 186 | ! |
---|
| 187 | CALL meter_to_levels(G1) |
---|
| 188 | ! |
---|
| 189 | DEALLOCATE(bathy_fin_constant) |
---|
| 190 | ! |
---|
| 191 | ! |
---|
| 192 | !------------------------------------------------------------------ |
---|
| 193 | ! ****Second option : new topo file or/and partial steps |
---|
| 194 | !------------------------------------------------------------------ |
---|
| 195 | ! |
---|
| 196 | ! |
---|
| 197 | ! |
---|
| 198 | ! |
---|
| 199 | ! |
---|
| 200 | ! |
---|
| 201 | ! |
---|
| 202 | ! |
---|
| 203 | ELSE |
---|
| 204 | ! |
---|
| 205 | WRITE(*,*) 'Second option : partial steps' |
---|
| 206 | ! read fine and coarse grids coordinates file |
---|
| 207 | ! |
---|
| 208 | status = Read_Coordinates(TRIM(parent_coordinate_file),G0) |
---|
| 209 | status = Read_Coordinates(TRIM(Child_coordinates),G1,Pacifique) |
---|
| 210 | ! |
---|
| 211 | IF( imax > SIZE(G0%nav_lon,1) .OR. jmax > SIZE(G0%nav_lon,2) .OR. & |
---|
| 212 | imax <= imin .OR. jmax <= jmin ) THEN |
---|
| 213 | WRITE(*,*) 'ERROR ***** bad child grid definition ...' |
---|
| 214 | WRITE(*,*) 'please check imin,jmin,imax,jmax,jpizoom,jpjzoom values' |
---|
| 215 | WRITE(*,*) ' ' |
---|
| 216 | STOP |
---|
| 217 | ENDIF |
---|
| 218 | ! |
---|
| 219 | |
---|
| 220 | ! |
---|
| 221 | IF( SIZE(G1%nav_lon,1) .NE. nxfin .OR. SIZE(G1%nav_lon,2) .NE. nyfin ) THEN |
---|
| 222 | ! |
---|
| 223 | WRITE(*,*) 'ERROR ***** bad child coordinates file ...' |
---|
| 224 | WRITE(*,*) ' ' |
---|
| 225 | WRITE(*,*) 'please check that child coordinates file ' |
---|
| 226 | WRITE(*,*) 'has been created with the same namelist ' |
---|
| 227 | WRITE(*,*) ' ' |
---|
| 228 | STOP |
---|
| 229 | ! |
---|
| 230 | ENDIF |
---|
| 231 | ! |
---|
| 232 | ! read coarse grid bathymetry file |
---|
| 233 | !--- |
---|
| 234 | IF( new_topo ) THEN |
---|
| 235 | WRITE(*,*) 'use new topo file : ',TRIM(elevation_database) |
---|
| 236 | DEALLOCATE( G0%nav_lon, G0%nav_lat ) |
---|
| 237 | status = Read_bathy_meter(TRIM(elevation_database),G0,G1,Pacifique) |
---|
| 238 | ELSE |
---|
| 239 | WRITE(*,*) 'no new topo file' |
---|
| 240 | status = Read_Bathymeter(TRIM(parent_bathy_meter),G0) |
---|
| 241 | IF(Pacifique) THEN |
---|
[9149] | 242 | WHERE(G0%nav_lon < 0.0001) |
---|
[2136] | 243 | G0%nav_lon = G0%nav_lon + 360. |
---|
| 244 | END WHERE |
---|
| 245 | ENDIF |
---|
| 246 | ENDIF |
---|
| 247 | !--- |
---|
| 248 | numlatfine = SIZE(G1%nav_lat,2) |
---|
| 249 | numlonfine = SIZE(G1%nav_lat,1) |
---|
| 250 | ! |
---|
| 251 | ! |
---|
| 252 | ALLOCATE(G1%bathy_meter(nxfin,nyfin)) |
---|
| 253 | G1%bathy_meter(:,:)=0. |
---|
| 254 | |
---|
| 255 | WRITE(*,*) 'Interpolation of high resolution bathymetry on child grid' |
---|
| 256 | |
---|
| 257 | IF( type_bathy_interp == 0 ) THEN |
---|
| 258 | WRITE(*,*) 'Arithmetic average ...' |
---|
| 259 | ELSE IF( type_bathy_interp == 1 ) THEN |
---|
| 260 | WRITE(*,*) 'Median average ...' |
---|
| 261 | ELSE IF( type_bathy_interp == 2 ) THEN |
---|
| 262 | WRITE(*,*) 'Bilinear interpolation ...' |
---|
| 263 | ELSE |
---|
| 264 | WRITE(*,*) 'bad value for type_bathy_interp variable ( must be 0,1 or 2 )' |
---|
| 265 | STOP |
---|
| 266 | ENDIF |
---|
| 267 | ! |
---|
| 268 | !************************************ |
---|
| 269 | !MEDIAN AVERAGE or ARITHMETIC AVERAGE |
---|
| 270 | !************************************ |
---|
| 271 | ! |
---|
| 272 | IF( type_bathy_interp == 0 .OR. type_bathy_interp == 1 ) THEN |
---|
| 273 | ! |
---|
| 274 | ALLOCATE(trouble_points(nxfin,nyfin)) |
---|
| 275 | trouble_points = 0 |
---|
| 276 | ! |
---|
| 277 | ! POINT DETECTION |
---|
| 278 | ! |
---|
| 279 | ! |
---|
| 280 | DO jj = 2,numlatfine |
---|
| 281 | DO ji = 2,numlonfine |
---|
| 282 | ! |
---|
| 283 | ! FINE GRID CELLS DEFINITION |
---|
| 284 | ! |
---|
| 285 | Cell_lonmin = MIN(G1%glamf(ji-1,jj-1),G1%glamf(ji,jj-1),G1%glamf(ji,jj),G1%glamf(ji-1,jj)) |
---|
| 286 | Cell_lonmax = MAX(G1%glamf(ji-1,jj-1),G1%glamf(ji,jj-1),G1%glamf(ji,jj),G1%glamf(ji-1,jj)) |
---|
| 287 | Cell_latmin = MIN(G1%gphif(ji-1,jj-1),G1%gphif(ji,jj-1),G1%gphif(ji,jj),G1%gphif(ji-1,jj)) |
---|
| 288 | Cell_latmax = MAX(G1%gphif(ji-1,jj-1),G1%gphif(ji,jj-1),G1%gphif(ji,jj),G1%gphif(ji-1,jj)) |
---|
| 289 | ! |
---|
| 290 | ! SEARCH FOR ALL POINTS CONTAINED IN THIS CELL |
---|
| 291 | ! |
---|
| 292 | iimin = 1 |
---|
| 293 | DO WHILE( G0%nav_lon(iimin,1) < Cell_lonmin ) |
---|
| 294 | iimin = iimin + 1 |
---|
| 295 | ENDDO |
---|
| 296 | ! |
---|
| 297 | jjmin = 1 |
---|
| 298 | DO WHILE( G0%nav_lat(iimin,jjmin) < Cell_latmin ) |
---|
| 299 | jjmin = jjmin + 1 |
---|
| 300 | ENDDO |
---|
| 301 | ! |
---|
| 302 | iimax = iimin |
---|
| 303 | DO WHILE( G0%nav_lon(iimax,1) <= Cell_lonmax ) |
---|
| 304 | iimax = iimax + 1 |
---|
[9149] | 305 | iimax = MIN( iimax,SIZE(G0%bathy_meter,1)) |
---|
[2136] | 306 | ENDDO |
---|
| 307 | ! |
---|
| 308 | jjmax = jjmin |
---|
| 309 | DO WHILE( G0%nav_lat(iimax,jjmax) <= Cell_latmax ) |
---|
| 310 | jjmax = jjmax + 1 |
---|
[9149] | 311 | jjmax = MIN( jjmax,SIZE(G0%bathy_meter,2)) |
---|
| 312 | |
---|
[2136] | 313 | ENDDO |
---|
| 314 | ! |
---|
| 315 | iimax = iimax-1 |
---|
| 316 | jjmax = jjmax-1 |
---|
| 317 | ! |
---|
| 318 | iimin = MAX( iimin,1 ) |
---|
| 319 | jjmin = MAX( jjmin,1 ) |
---|
| 320 | iimax = MIN( iimax,SIZE(G0%bathy_meter,1)) |
---|
| 321 | jjmax = MIN( jjmax,SIZE(G0%bathy_meter,2)) |
---|
| 322 | |
---|
| 323 | nxhr = iimax - iimin + 1 |
---|
| 324 | nyhr = jjmax - jjmin + 1 |
---|
| 325 | |
---|
| 326 | IF( nxhr == 0 .OR. nyhr == 0 ) THEN |
---|
| 327 | trouble_points(ji,jj) = 1 |
---|
| 328 | ELSE |
---|
| 329 | |
---|
| 330 | ALLOCATE( vardep(nxhr,nyhr) ) |
---|
| 331 | ALLOCATE( mask_oce(nxhr,nyhr) ) |
---|
| 332 | mask_oce = 0 |
---|
| 333 | |
---|
| 334 | vardep(:,:) = G0%bathy_meter(iimin:iimax,jjmin:jjmax) |
---|
| 335 | |
---|
| 336 | WHERE( vardep(:,:) .GT. 0. ) mask_oce = 1 |
---|
| 337 | |
---|
| 338 | IF( SUM(mask_oce) == 0 ) THEN |
---|
| 339 | G1%bathy_meter(ji,jj) = 0. |
---|
| 340 | ELSE |
---|
| 341 | IF( type_bathy_interp == 0 ) THEN |
---|
| 342 | ! Arithmetic average |
---|
| 343 | G1%bathy_meter(ji,jj) = SUM (vardep(:,:)*mask_oce(:,:))/SUM(mask_oce) |
---|
| 344 | ELSE |
---|
| 345 | ! Median average |
---|
| 346 | ! |
---|
| 347 | vardep(:,:) = vardep(:,:)*mask_oce(:,:) - 100000*(1-mask_oce(:,:)) |
---|
| 348 | ALLOCATE(vardep1d(nxhr*nyhr)) |
---|
| 349 | vardep1d = RESHAPE(vardep,(/ nxhr*nyhr /) ) |
---|
| 350 | CALL ssort(vardep1d,nxhr*nyhr) |
---|
| 351 | ! |
---|
| 352 | ! Calculate median |
---|
| 353 | ! |
---|
| 354 | IF (MOD(SUM(mask_oce),2) .NE. 0) THEN |
---|
| 355 | G1%bathy_meter(ji,jj) = vardep1d( SUM(mask_oce)/2 + 1) |
---|
| 356 | ELSE |
---|
| 357 | G1%bathy_meter(ji,jj) = ( vardep1d(SUM(mask_oce)/2) + vardep1d(SUM(mask_oce)/2+1) )/2.0 |
---|
| 358 | END IF |
---|
| 359 | DEALLOCATE(vardep1d) |
---|
| 360 | ENDIF |
---|
| 361 | ENDIF |
---|
| 362 | DEALLOCATE (mask_oce,vardep) |
---|
| 363 | |
---|
| 364 | ENDIF |
---|
| 365 | ENDDO |
---|
| 366 | ENDDO |
---|
| 367 | |
---|
| 368 | IF( SUM( trouble_points ) > 0 ) THEN |
---|
| 369 | PRINT*,'too much empty cells, proceed to bilinear interpolation !!' |
---|
| 370 | type_bathy_interp = 2 |
---|
| 371 | ENDIF |
---|
| 372 | |
---|
| 373 | ENDIF |
---|
| 374 | |
---|
| 375 | ! |
---|
| 376 | ! create logical array masksrc |
---|
| 377 | ! |
---|
| 378 | IF( type_bathy_interp == 2) THEN |
---|
| 379 | ! |
---|
| 380 | |
---|
| 381 | ! |
---|
| 382 | identical_grids = .FALSE. |
---|
| 383 | |
---|
| 384 | IF( SIZE(G0%nav_lat,1) == SIZE(G1%nav_lat,1) .AND. & |
---|
| 385 | SIZE(G0%nav_lat,2) == SIZE(G1%nav_lat,2) .AND. & |
---|
| 386 | SIZE(G0%nav_lon,1) == SIZE(G1%nav_lon,1) .AND. & |
---|
| 387 | SIZE(G0%nav_lon,2) == SIZE(G1%nav_lon,2) ) THEN |
---|
| 388 | IF( MAXVAL( ABS(G0%nav_lat(:,:)- G1%nav_lat(:,:)) ) < 0.0001 .AND. & |
---|
| 389 | MAXVAL( ABS(G0%nav_lon(:,:)- G1%nav_lon(:,:)) ) < 0.0001 ) THEN |
---|
| 390 | G1%bathy_meter = G0%bathy_meter |
---|
| 391 | PRINT*,'same grid between ',elevation_database,' and child domain' |
---|
| 392 | identical_grids = .TRUE. |
---|
| 393 | ENDIF |
---|
| 394 | ENDIF |
---|
| 395 | |
---|
| 396 | |
---|
| 397 | IF( .NOT. identical_grids ) THEN |
---|
| 398 | |
---|
| 399 | ALLOCATE(masksrc(SIZE(G0%bathy_meter,1),SIZE(G0%bathy_meter,2))) |
---|
| 400 | ALLOCATE(bathy_test(nxfin,nyfin)) |
---|
| 401 | ! |
---|
| 402 | ! Where(G0%bathy_meter.le.0.00001) |
---|
| 403 | ! masksrc = .false. |
---|
| 404 | ! ElseWhere |
---|
| 405 | ! |
---|
| 406 | masksrc = .TRUE. |
---|
| 407 | ! |
---|
| 408 | ! End where |
---|
| 409 | ! |
---|
| 410 | ! compute remapping matrix thanks to SCRIP package |
---|
| 411 | ! |
---|
| 412 | CALL get_remap_matrix(G0%nav_lat,G1%nav_lat, & |
---|
| 413 | G0%nav_lon,G1%nav_lon, & |
---|
| 414 | masksrc,matrix,src_add,dst_add) |
---|
| 415 | CALL make_remap(G0%bathy_meter,bathy_test,nxfin,nyfin, & |
---|
| 416 | matrix,src_add,dst_add) |
---|
| 417 | ! |
---|
| 418 | G1%bathy_meter = bathy_test |
---|
| 419 | ! |
---|
| 420 | DEALLOCATE(masksrc) |
---|
| 421 | DEALLOCATE(bathy_test) |
---|
| 422 | |
---|
| 423 | ENDIF |
---|
| 424 | ! |
---|
| 425 | ENDIF |
---|
| 426 | ! |
---|
| 427 | ! |
---|
| 428 | ! |
---|
| 429 | !------------------------------------------------------------------------------------------ |
---|
| 430 | ! ! ****Third option : Partial Steps |
---|
| 431 | ! The code includes the |
---|
| 432 | ! option to include partial cells at the bottom |
---|
| 433 | ! in order to better resolve topographic variations |
---|
| 434 | !------------------------------------------------------------------------------------------ |
---|
| 435 | ! |
---|
| 436 | ! At this step bathymetry in meters has already been interpolated on fine grid |
---|
| 437 | ! |
---|
| 438 | ! |
---|
| 439 | IF( partial_steps ) THEN |
---|
| 440 | ! |
---|
| 441 | status = Read_Bathymeter(TRIM(parent_bathy_meter),G0) |
---|
| 442 | DEALLOCATE(G0%nav_lat,G0%nav_lon) |
---|
| 443 | status = Read_coordinates(TRIM(parent_coordinate_file),G0) |
---|
| 444 | !------------------------ |
---|
| 445 | IF (.NOT.ASSOCIATED(G0%gdepw_ps)) & |
---|
| 446 | ALLOCATE(G0%gdepw_ps(SIZE(G0%bathy_meter,1),SIZE(G0%bathy_meter,2))) |
---|
| 447 | IF (.NOT.ASSOCIATED(G1%gdepw_ps)) & |
---|
| 448 | ALLOCATE(G1%gdepw_ps(SIZE(G1%bathy_meter,1),SIZE(G1%bathy_meter,2))) |
---|
| 449 | IF (.NOT.ASSOCIATED(gdepw_ps_interp)) & |
---|
| 450 | ALLOCATE(gdepw_ps_interp(SIZE(G1%bathy_meter,1),SIZE(G1%bathy_meter,2))) |
---|
| 451 | ! |
---|
| 452 | ! |
---|
| 453 | WRITE(*,*) 'Coarse grid : ' |
---|
| 454 | CALL get_partial_steps(G0) |
---|
| 455 | WRITE(*,*) ' ' |
---|
| 456 | WRITE(*,*) 'Fine grid : ' |
---|
| 457 | CALL get_partial_steps(G1) ! compute gdepw_ps for G1 |
---|
| 458 | CALL bathymetry_control(G0%Bathy_level) ! |
---|
| 459 | CALL Check_interp(G0,gdepw_ps_interp) ! interpolation in connection zone (3 coarse grid cells) |
---|
| 460 | ! |
---|
[9166] | 461 | boundary = connectionsize*irafx + nbghostcellsfine + 1 |
---|
[2136] | 462 | G1%gdepw_ps(1:boundary,:) = gdepw_ps_interp(1:boundary,:) |
---|
| 463 | G1%gdepw_ps(:,1:boundary) = gdepw_ps_interp(:,1:boundary) |
---|
| 464 | G1%gdepw_ps(nxfin-boundary+1:nxfin,:) = gdepw_ps_interp(nxfin-boundary+1:nxfin,:) |
---|
| 465 | G1%gdepw_ps(:,nyfin-boundary+1:nyfin) = gdepw_ps_interp(:,nyfin-boundary+1:nyfin) |
---|
| 466 | |
---|
| 467 | |
---|
| 468 | ! |
---|
| 469 | |
---|
| 470 | IF(.NOT. smoothing) WRITE(*,*) 'No smoothing process only connection is carried out' |
---|
| 471 | WRITE(*,*) ' linear connection on ',nb_connection_pts,'coarse grid points' |
---|
| 472 | |
---|
[9149] | 473 | connectionsize = 2 + nb_connection_pts |
---|
[2136] | 474 | ! |
---|
| 475 | gdepw_ps_interp = 0. |
---|
| 476 | CALL Check_interp(G0,gdepw_ps_interp) ! interpolation in connection zone (3 coarse grid cells) |
---|
| 477 | ! |
---|
| 478 | ! |
---|
| 479 | ! |
---|
| 480 | ! |
---|
| 481 | ! LINEAR CONNECTION |
---|
| 482 | ! |
---|
| 483 | ! |
---|
| 484 | ! |
---|
| 485 | ! |
---|
| 486 | ! |
---|
[9149] | 487 | IF (.NOT.ASSOCIATED(wtab)) & |
---|
| 488 | ALLOCATE(wtab(SIZE(G1%bathy_meter,1),SIZE(G1%bathy_meter,2))) |
---|
| 489 | wtab(:,:) = 0. |
---|
[2136] | 490 | wghts = 1. |
---|
| 491 | DO ji = boundary + 1 , boundary + nb_connection_pts * irafx |
---|
[9149] | 492 | wghts = wghts - (1. / (nb_connection_pts*irafx - 1. ) ) |
---|
| 493 | DO jj=boundary+1,nyfin-boundary |
---|
| 494 | wtab(ji,jj) = MAX(wghts, wtab(ji,jj)) |
---|
| 495 | END DO |
---|
| 496 | END DO |
---|
| 497 | |
---|
| 498 | wghts = 1. |
---|
| 499 | DO ji = nxfin - boundary , nxfin - boundary - nb_connection_pts * irafx +1 , -1 |
---|
| 500 | wghts = wghts - (1. / (nb_connection_pts*irafx - 1. ) ) |
---|
| 501 | DO jj=boundary+1,nyfin-boundary |
---|
| 502 | wtab(ji,jj) = MAX(wghts, wtab(ji,jj)) |
---|
| 503 | END DO |
---|
| 504 | END DO |
---|
[2136] | 505 | |
---|
| 506 | wghts = 1. |
---|
| 507 | DO jj = boundary + 1 , boundary + nb_connection_pts * irafy |
---|
[9149] | 508 | wghts = wghts - (1. / (nb_connection_pts*irafy - 1. ) ) |
---|
| 509 | DO ji=boundary+1,nxfin-boundary |
---|
| 510 | wtab(ji,jj) = MAX(wghts, wtab(ji,jj)) |
---|
| 511 | END DO |
---|
| 512 | END DO |
---|
| 513 | |
---|
[2136] | 514 | wghts = 1. |
---|
[9149] | 515 | DO jj = nyfin - boundary , nyfin - boundary - nb_connection_pts * irafy +1, -1 |
---|
| 516 | wghts = wghts - (1. / (nb_connection_pts*irafy - 1. ) ) |
---|
| 517 | DO ji=boundary+1,nxfin-boundary |
---|
| 518 | wtab(ji,jj) = MAX(wghts, wtab(ji,jj)) |
---|
| 519 | END DO |
---|
| 520 | END DO |
---|
[2136] | 521 | |
---|
[9149] | 522 | G1%gdepw_ps(:,:) = (1.-wtab(:,:)) * G1%gdepw_ps(:,:) + gdepw_ps_interp(:,:)*wtab(:,:) |
---|
| 523 | |
---|
[2136] | 524 | G1%bathy_meter = G1%gdepw_ps |
---|
| 525 | ! |
---|
[9149] | 526 | connectionsize = 2 |
---|
[2136] | 527 | ! |
---|
| 528 | IF(smoothing) THEN |
---|
| 529 | |
---|
| 530 | ! |
---|
| 531 | ! Smoothing to connect the connection zone (3 + nb_connection_pts coarse grid cells) and the interior domain |
---|
| 532 | ! |
---|
[9166] | 533 | boundary = (connectionsize+nb_connection_pts)*irafx + nbghostcellsfine + 1 |
---|
[2136] | 534 | CALL smooth_topo(G1%gdepw_ps(boundary:nxfin-boundary+1,boundary:nyfin-boundary+1),nbiter) |
---|
| 535 | G1%bathy_meter = G1%gdepw_ps |
---|
| 536 | ENDIF |
---|
| 537 | |
---|
| 538 | |
---|
| 539 | ! |
---|
| 540 | ! |
---|
| 541 | ! Remove closed seas |
---|
| 542 | ! |
---|
| 543 | IF (removeclosedseas) THEN |
---|
| 544 | ALLOCATE(bathy_test(nxfin,nyfin)) |
---|
| 545 | bathy_test=0. |
---|
| 546 | WHERE (G1%bathy_meter(1,:).GT.0.) |
---|
| 547 | bathy_test(1,:)=1 |
---|
| 548 | END WHERE |
---|
| 549 | WHERE (G1%bathy_meter(nxfin,:).GT.0.) |
---|
| 550 | bathy_test(nxfin,:)=1 |
---|
| 551 | END WHERE |
---|
| 552 | WHERE (G1%bathy_meter(:,1).GT.0.) |
---|
| 553 | bathy_test(:,1)=1 |
---|
| 554 | END WHERE |
---|
| 555 | WHERE (G1%bathy_meter(:,nyfin).GT.0.) |
---|
| 556 | bathy_test(:,nyfin)=1 |
---|
| 557 | END WHERE |
---|
| 558 | nbadd = 1 |
---|
| 559 | DO WHILE (nbadd.NE.0) |
---|
| 560 | nbadd = 0 |
---|
| 561 | DO j=2,nyfin-1 |
---|
| 562 | DO i=2,nxfin-1 |
---|
| 563 | IF (G1%bathy_meter(i,j).GT.0.) THEN |
---|
| 564 | IF (MAX(bathy_test(i,j+1),bathy_test(i,j-1), & |
---|
| 565 | bathy_test(i-1,j),bathy_test(i+1,j)).EQ.1) THEN |
---|
| 566 | IF (bathy_test(i,j).NE.1.) nbadd = nbadd + 1 |
---|
| 567 | bathy_test(i,j)=1. |
---|
| 568 | ENDIF |
---|
| 569 | |
---|
| 570 | ENDIF |
---|
| 571 | ENDDO |
---|
| 572 | ENDDO |
---|
| 573 | ENDDO |
---|
| 574 | WHERE (bathy_test.EQ.0.) |
---|
| 575 | G1%bathy_meter = 0. |
---|
| 576 | END WHERE |
---|
| 577 | DEALLOCATE(bathy_test) |
---|
| 578 | ENDIF |
---|
| 579 | ! |
---|
| 580 | IF(bathy_update) CALL Update_Parent_Bathy( G0,G1 ) |
---|
| 581 | ! |
---|
| 582 | CALL set_child_name(parent_bathy_meter,child_ps) |
---|
| 583 | status = Write_Bathy_meter(TRIM(child_ps),G1) |
---|
| 584 | |
---|
| 585 | IF(bathy_update) status = Write_Bathy_meter(TRIM(updated_parent_file),G0) |
---|
| 586 | |
---|
| 587 | CALL get_partial_steps(G1) |
---|
| 588 | ! |
---|
| 589 | G1%bathy_level=NINT(G1%bathy_level) |
---|
| 590 | ! |
---|
[2455] | 591 | IF( TRIM(parent_meshmask_file) .NE. '/NULL' ) & |
---|
[2136] | 592 | status = Write_Bathy_level(TRIM(Childlevel_file),G1) |
---|
| 593 | ! |
---|
| 594 | WRITE(*,*) '****** Bathymetry successfully created for partial cells ******' |
---|
| 595 | WRITE(*,*) ' ' |
---|
| 596 | ! |
---|
| 597 | STOP |
---|
| 598 | ENDIF |
---|
| 599 | ! |
---|
| 600 | !--------------------------------end if partial steps------------------------------------ |
---|
| 601 | ! |
---|
| 602 | ! |
---|
[2455] | 603 | status = Read_bathy_level(TRIM(parent_meshmask_file),G0) |
---|
[2136] | 604 | ! |
---|
| 605 | CALL levels_to_meter(G0) |
---|
| 606 | ! |
---|
| 607 | ! compute constant bathymetry for Parent-Child bathymetry connection |
---|
| 608 | ! |
---|
| 609 | WHERE( G0%bathy_meter < 0.000001 ) G0%bathy_meter = 0. |
---|
| 610 | |
---|
| 611 | CALL init_constant_bathy(G0%bathy_meter,bathy_fin_constant) |
---|
| 612 | ! |
---|
[9166] | 613 | boundary = connectionsize*irafx + nbghostcellsfine + 1 |
---|
[2136] | 614 | ! |
---|
| 615 | G1%bathy_meter(1:boundary,:) = bathy_fin_constant(1:boundary,:) |
---|
| 616 | G1%bathy_meter(:,1:boundary) = bathy_fin_constant(:,1:boundary) |
---|
| 617 | G1%bathy_meter(nxfin-boundary+1:nxfin,:) = bathy_fin_constant(nxfin-boundary+1:nxfin,:) |
---|
| 618 | G1%bathy_meter(:,nyfin-boundary+1:nyfin) = bathy_fin_constant(:,nyfin-boundary+1:nyfin) |
---|
| 619 | ! |
---|
| 620 | ! bathymetry smoothing |
---|
| 621 | ! |
---|
| 622 | CALL smooth_topo(G1%bathy_meter(boundary:nxfin-boundary+1,boundary:nyfin-boundary+1),nbiter) |
---|
| 623 | ! |
---|
| 624 | ! convert bathymetry from meters to levels |
---|
| 625 | ! |
---|
| 626 | CALL meter_to_levels(G1) |
---|
| 627 | ! |
---|
| 628 | DEALLOCATE(G1%bathy_meter) |
---|
| 629 | ! |
---|
| 630 | ENDIF |
---|
| 631 | ! |
---|
| 632 | ! |
---|
| 633 | ! make connection thanks to constant and interpolated bathymetry |
---|
| 634 | ! |
---|
| 635 | ! |
---|
| 636 | G1%bathy_level=NINT(G1%bathy_level) |
---|
| 637 | ! |
---|
| 638 | DO j=1,nyfin |
---|
| 639 | DO i=1,nxfin |
---|
| 640 | IF (g1%bathy_level(i,j).LT.0.) THEN |
---|
| 641 | PRINT *,'error in ',i,j,g1%bathy_level(i,j) |
---|
| 642 | ENDIF |
---|
| 643 | ENDDO |
---|
| 644 | ENDDO |
---|
| 645 | ! |
---|
| 646 | WHERE ((G1%bathy_level.LT.3.).AND.(G1%bathy_level.GT.0.)) |
---|
| 647 | G1%bathy_level=3 |
---|
| 648 | END WHERE |
---|
| 649 | ! |
---|
| 650 | ! possibility to remove closed seas |
---|
| 651 | ! |
---|
| 652 | IF (removeclosedseas) THEN |
---|
| 653 | ALLOCATE(bathy_test(nxfin,nyfin)) |
---|
| 654 | |
---|
| 655 | bathy_test=0. |
---|
| 656 | WHERE (G1%bathy_level(1,:).GT.0.) |
---|
| 657 | bathy_test(1,:)=1 |
---|
| 658 | END WHERE |
---|
| 659 | |
---|
| 660 | WHERE (G1%bathy_level(nxfin,:).GT.0.) |
---|
| 661 | bathy_test(nxfin,:)=1 |
---|
| 662 | END WHERE |
---|
| 663 | |
---|
| 664 | WHERE (G1%bathy_level(:,1).GT.0.) |
---|
| 665 | bathy_test(:,1)=1 |
---|
| 666 | END WHERE |
---|
| 667 | |
---|
| 668 | WHERE (G1%bathy_level(:,nyfin).GT.0.) |
---|
| 669 | bathy_test(:,nyfin)=1 |
---|
| 670 | END WHERE |
---|
| 671 | |
---|
| 672 | nbadd = 1 |
---|
| 673 | |
---|
| 674 | DO WHILE (nbadd.NE.0) |
---|
| 675 | nbadd = 0 |
---|
| 676 | DO j=2,nyfin-1 |
---|
| 677 | DO i=2,nxfin-1 |
---|
| 678 | IF (G1%bathy_level(i,j).GT.0.) THEN |
---|
| 679 | IF (MAX(bathy_test(i,j+1),bathy_test(i,j-1),bathy_test(i-1,j),bathy_test(i+1,j)).EQ.1) THEN |
---|
| 680 | IF (bathy_test(i,j).NE.1.) nbadd = nbadd + 1 |
---|
| 681 | bathy_test(i,j)=1. |
---|
| 682 | ENDIF |
---|
| 683 | |
---|
| 684 | ENDIF |
---|
| 685 | ENDDO |
---|
| 686 | ENDDO |
---|
| 687 | |
---|
| 688 | ENDDO |
---|
| 689 | |
---|
| 690 | WHERE (bathy_test.EQ.0.) |
---|
| 691 | G1%bathy_level = 0. |
---|
| 692 | END WHERE |
---|
| 693 | DEALLOCATE(bathy_test) |
---|
| 694 | ENDIF |
---|
| 695 | |
---|
| 696 | |
---|
| 697 | ! |
---|
| 698 | ! store interpolation result in output file |
---|
| 699 | ! |
---|
| 700 | status = Write_Bathy_level(TRIM(Childlevel_file),G1) |
---|
| 701 | |
---|
| 702 | WRITE(*,*) '****** Bathymetry successfully created for full cells ******' |
---|
| 703 | WRITE(*,*) ' ' |
---|
| 704 | |
---|
| 705 | STOP |
---|
| 706 | END PROGRAM create_bathy |
---|
| 707 | |
---|
| 708 | |
---|