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