[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 |
---|
[10025] | 15 | USE agrif_interpolation |
---|
[2136] | 16 | ! |
---|
| 17 | IMPLICIT NONE |
---|
| 18 | ! |
---|
| 19 | !************************************************************************ |
---|
| 20 | ! * |
---|
[2455] | 21 | ! PROGRAM CREATE_BATHY * |
---|
[2136] | 22 | ! * |
---|
| 23 | ! program to implement bathymetry interpolation to generate * |
---|
| 24 | ! child grid bathymetry file * |
---|
| 25 | ! * |
---|
| 26 | ! various options : * |
---|
| 27 | ! * |
---|
| 28 | ! 1- Interpolation directly from parent bathymetry file (z-coord) * |
---|
[2455] | 29 | ! 2- Use new topo file in meters (for example etopo) * |
---|
[2136] | 30 | ! * |
---|
| 31 | ! vertical coordinates permitted : z-coord and partial steps * |
---|
| 32 | ! sigma coordinates is not yet implemented * |
---|
| 33 | ! * |
---|
| 34 | !Interpolation is carried out using bilinear interpolation * |
---|
| 35 | !routine from SCRIP package or median average * |
---|
| 36 | ! * |
---|
| 37 | !http://climate.lanl.gov/Software/SCRIP/ * |
---|
| 38 | !************************************************************************ |
---|
| 39 | ! |
---|
| 40 | ! variables declaration |
---|
| 41 | ! |
---|
[10027] | 42 | CHARACTER(len=80) :: namelistname |
---|
[10381] | 43 | CHARACTER*100 :: child_coordinates, child_level, child_meter, child_domcfg |
---|
[10027] | 44 | LOGICAL :: identical_grids |
---|
| 45 | INTEGER :: nbadd,status,narg,iargc |
---|
| 46 | INTEGER :: jpj,jpi |
---|
| 47 | LOGICAL,DIMENSION(:,:),POINTER :: masksrc => NULL() |
---|
| 48 | INTEGER,DIMENSION(:,:),ALLOCATABLE :: mask_oce,trouble_points |
---|
| 49 | INTEGER,DIMENSION(:) ,POINTER :: src_add,dst_add => NULL() |
---|
| 50 | REAL*8, DIMENSION(:,:),POINTER :: matrix,interpdata => NULL() |
---|
| 51 | REAL*8, DIMENSION(:,:),POINTER :: bathy_fin_constant => NULL() |
---|
| 52 | REAL*8, DIMENSION(:,:),ALLOCATABLE :: bathy_test,vardep |
---|
| 53 | REAL*8, DIMENSION(:) ,ALLOCATABLE :: vardep1d |
---|
| 54 | REAL*8, DIMENSION(:,:),POINTER :: gdepw_ps_interp => NULL() |
---|
| 55 | REAL*8 :: Cell_lonmin,Cell_lonmax,Cell_latmin,Cell_latmax,wghts |
---|
| 56 | LOGICAL :: Pacifique = .FALSE. |
---|
| 57 | INTEGER :: boundary,iimin,iimax,jjmax,jjmin |
---|
| 58 | INTEGER :: nxhr,nyhr,nxyhr,ji,jj,nbiter |
---|
[2136] | 59 | |
---|
| 60 | TYPE(Coordinates) :: G0,G1 |
---|
| 61 | ! |
---|
| 62 | narg = iargc() |
---|
| 63 | IF (narg == 0) THEN |
---|
| 64 | namelistname = 'namelist.input' |
---|
| 65 | ELSE |
---|
| 66 | CALL getarg(1,namelistname) |
---|
| 67 | ENDIF |
---|
| 68 | ! |
---|
| 69 | ! read input file (namelist.input) |
---|
| 70 | CALL read_namelist(namelistname) |
---|
[10381] | 71 | |
---|
| 72 | ! if level or meter name is missing |
---|
| 73 | IF( TRIM(parent_level_name) == '' ) parent_level_name='mbathy' |
---|
| 74 | IF( TRIM(parent_meter_name) == '' ) parent_meter_name='Bathymetry' |
---|
| 75 | |
---|
[2136] | 76 | ! define names of child grid files |
---|
[10381] | 77 | CALL set_child_name(parent_coordinate_file,child_coordinates) |
---|
| 78 | IF( TRIM(parent_bathy_level) /= '' ) CALL set_child_name(parent_bathy_level,child_level) |
---|
| 79 | IF( TRIM(parent_bathy_meter) /= '' ) CALL set_child_name(parent_bathy_meter,child_meter) |
---|
| 80 | IF( TRIM(parent_domcfg_out) /= '' ) CALL set_child_name(parent_domcfg_out,child_domcfg) |
---|
[2136] | 81 | ! |
---|
[10381] | 82 | IF( TRIM(parent_bathy_level) == '' .AND. TRIM(parent_bathy_meter) == '') THEN |
---|
| 83 | WRITE(*,*) 'ERROR ***** one needs at least to define parent_bathy_level or parent_bathy_meter ...' |
---|
| 84 | STOP |
---|
| 85 | ENDIF |
---|
[10025] | 86 | ! ! ------------------------------------------------------------------ |
---|
| 87 | IF( .NOT.new_topo .AND. .NOT.partial_steps ) THEN ! **** First option : no new topo file & no partial steps |
---|
| 88 | ! ! ==> interpolate levels directly from parent file |
---|
| 89 | ! ! ------------------------------------------------------------------ |
---|
| 90 | WRITE(*,*) '*** First option: no new topo file & no partial steps' |
---|
[2136] | 91 | ! |
---|
[10381] | 92 | ! read coarse grid coordinates |
---|
[2136] | 93 | status = Read_Coordinates(TRIM(parent_coordinate_file),G0) |
---|
[10381] | 94 | |
---|
| 95 | ! read coarse grid bathymetry |
---|
| 96 | IF( TRIM(parent_bathy_level) /= '' ) THEN |
---|
| 97 | status = Read_bathy_level(TRIM(parent_bathy_level),G0) |
---|
| 98 | CALL levels_to_meter(G0) |
---|
| 99 | ELSEIF( TRIM(parent_bathy_level) == '' .AND. TRIM(parent_bathy_meter) /= '') THEN |
---|
| 100 | status = read_bathy_meter(TRIM(parent_bathy_meter),G0) |
---|
| 101 | CALL meter_to_levels(G0) |
---|
| 102 | ENDIF |
---|
[2136] | 103 | ! |
---|
[10025] | 104 | ! read fine grid coordinates |
---|
[10381] | 105 | status = Read_Coordinates(TRIM(child_coordinates),G1,pacifique) |
---|
[10025] | 106 | ! |
---|
| 107 | ! stop if error in size |
---|
| 108 | IF( imax > SIZE(G0%glamt,1) .OR. jmax > SIZE(G0%glamt,2) .OR. imax <= imin .OR. jmax <= jmin ) THEN |
---|
[2136] | 109 | WRITE(*,*) 'ERROR ***** bad child grid definition ...' |
---|
| 110 | WRITE(*,*) 'please check imin,jmin,imax,jmax,jpizoom,jpjzoom values' |
---|
| 111 | STOP |
---|
| 112 | ENDIF |
---|
| 113 | IF( SIZE(G1%nav_lon,1) .NE. nxfin .OR. SIZE(G1%nav_lon,2) .NE. nyfin ) THEN |
---|
| 114 | WRITE(*,*) 'ERROR ***** bad child coordinates file ...' |
---|
[10025] | 115 | WRITE(*,*) 'please check that child coordinates file has been created with the same namelist' |
---|
[2136] | 116 | STOP |
---|
| 117 | ENDIF |
---|
| 118 | ! |
---|
[10025] | 119 | jpj = SIZE(G0%nav_lat,2) |
---|
| 120 | jpi = SIZE(G0%nav_lat,1) |
---|
[2136] | 121 | ! |
---|
| 122 | ! create logical array masksrc |
---|
[10025] | 123 | ALLOCATE(masksrc(jpi,jpj)) |
---|
| 124 | WHERE(G0%bathy_level.LE.0) ; masksrc = .FALSE. ; |
---|
| 125 | ELSEWHERE ; masksrc = .TRUE. ; |
---|
[2136] | 126 | END WHERE |
---|
| 127 | |
---|
[10025] | 128 | ! change longitudes (from -180:180 to 0:360) |
---|
[2136] | 129 | IF ( Pacifique ) THEN |
---|
[10025] | 130 | WHERE(G0%nav_lon < 0.001) G0%nav_lon = G0%nav_lon + 360. |
---|
[2136] | 131 | ENDIF |
---|
[10381] | 132 | ! |
---|
[2136] | 133 | ALLOCATE(G1%bathy_meter(nxfin,nyfin)) |
---|
[10381] | 134 | ! |
---|
[10025] | 135 | ! compute remapping matrix thanks to SCRIP package |
---|
| 136 | CALL get_remap_matrix(G0%nav_lat,G1%nav_lat,G0%nav_lon,G1%nav_lon,masksrc,matrix,src_add,dst_add) |
---|
| 137 | CALL make_remap(G0%bathy_meter,G1%bathy_meter,nxfin,nyfin,matrix,src_add,dst_add) |
---|
[2136] | 138 | DEALLOCATE(masksrc) |
---|
| 139 | ! |
---|
| 140 | ! compute constant bathymetry for Parent-Child bathymetry connection |
---|
| 141 | CALL init_constant_bathy(G0%bathy_meter,bathy_fin_constant) |
---|
| 142 | ! |
---|
[10025] | 143 | ! replace child bathymetry by parent bathymetry at the boundaries |
---|
| 144 | IF( ln_agrif_domain ) THEN |
---|
| 145 | boundary = npt_copy*irafx + nbghostcellsfine + 1 |
---|
| 146 | ELSE |
---|
| 147 | boundary = npt_copy*irafx |
---|
| 148 | ENDIF |
---|
[2136] | 149 | G1%bathy_meter(1:boundary,:) = bathy_fin_constant(1:boundary,:) |
---|
| 150 | G1%bathy_meter(:,1:boundary) = bathy_fin_constant(:,1:boundary) |
---|
| 151 | G1%bathy_meter(nxfin-boundary+1:nxfin,:) = bathy_fin_constant(nxfin-boundary+1:nxfin,:) |
---|
| 152 | G1%bathy_meter(:,nyfin-boundary+1:nyfin) = bathy_fin_constant(:,nyfin-boundary+1:nyfin) |
---|
[10025] | 153 | DEALLOCATE(bathy_fin_constant) |
---|
[2136] | 154 | ! |
---|
[10025] | 155 | ! bathymetry smoothing (everywhere except at the boundaries) |
---|
| 156 | IF( smoothing ) THEN |
---|
| 157 | IF( ln_agrif_domain ) THEN |
---|
| 158 | CALL smooth_topo(G1%bathy_meter(boundary:nxfin-boundary+1,boundary:nyfin-boundary+1),nbiter) |
---|
| 159 | ELSE |
---|
| 160 | CALL smooth_topo(G1%bathy_meter(boundary+1:nxfin-boundary,boundary+1:nyfin-boundary),nbiter) |
---|
| 161 | ENDIF |
---|
| 162 | ENDIF |
---|
| 163 | ! |
---|
| 164 | ! From meters to levels |
---|
[2136] | 165 | CALL meter_to_levels(G1) |
---|
[10025] | 166 | G1%bathy_level=NINT(G1%bathy_level) |
---|
| 167 | ! |
---|
| 168 | ! Check errors in bathy |
---|
[10027] | 169 | DO jj=1,nyfin |
---|
| 170 | DO ji=1,nxfin |
---|
| 171 | IF (G1%bathy_level(ji,jj).LT.0.) THEN |
---|
| 172 | PRINT *,'error in ',ji,jj,G1%bathy_level(ji,jj) |
---|
[10025] | 173 | STOP |
---|
| 174 | ENDIF |
---|
| 175 | ENDDO |
---|
| 176 | ENDDO |
---|
| 177 | WHERE ((G1%bathy_level.LT.3.).AND.(G1%bathy_level.GT.0.)) G1%bathy_level=3 |
---|
[2136] | 178 | ! |
---|
[10025] | 179 | ! remove closed seas |
---|
| 180 | IF (removeclosedseas) THEN |
---|
| 181 | ALLOCATE(bathy_test(nxfin,nyfin)) |
---|
| 182 | |
---|
| 183 | bathy_test=0. |
---|
| 184 | WHERE (G1%bathy_level(1,:) .GT.0.) bathy_test(1,:)=1 |
---|
| 185 | WHERE (G1%bathy_level(nxfin,:) .GT.0.) bathy_test(nxfin,:)=1 |
---|
| 186 | WHERE (G1%bathy_level(:,1) .GT.0.) bathy_test(:,1)=1 |
---|
| 187 | WHERE (G1%bathy_level(:,nyfin) .GT.0.) bathy_test(:,nyfin)=1 |
---|
| 188 | |
---|
| 189 | nbadd = 1 |
---|
| 190 | DO WHILE (nbadd.NE.0) |
---|
| 191 | nbadd = 0 |
---|
[10027] | 192 | DO jj=2,nyfin-1 |
---|
| 193 | DO ji=2,nxfin-1 |
---|
| 194 | IF (G1%bathy_level(ji,jj).GT.0.) THEN |
---|
| 195 | IF (MAX(bathy_test(ji,jj+1),bathy_test(ji,jj-1),bathy_test(ji-1,jj),bathy_test(ji+1,jj)).EQ.1) THEN |
---|
| 196 | IF (bathy_test(ji,jj).NE.1.) nbadd = nbadd + 1 |
---|
| 197 | bathy_test(ji,jj)=1. |
---|
[10025] | 198 | ENDIF |
---|
| 199 | ENDIF |
---|
| 200 | ENDDO |
---|
| 201 | ENDDO |
---|
| 202 | ENDDO |
---|
| 203 | |
---|
| 204 | WHERE (bathy_test.EQ.0.) G1%bathy_level = 0. |
---|
| 205 | |
---|
[10381] | 206 | DEALLOCATE(bathy_test) |
---|
[10025] | 207 | ENDIF |
---|
[2136] | 208 | ! |
---|
[10381] | 209 | CALL levels_to_meter(G1) ! needed for domcfg |
---|
| 210 | ! |
---|
[10025] | 211 | ! store interpolation result in output file |
---|
[10381] | 212 | status = Write_Bathy_level(TRIM(child_level),G1) |
---|
| 213 | status = write_domcfg(TRIM(child_domcfg),G1) |
---|
| 214 | !!status = write_domcfg(TRIM(parent_domcfg_out),G0) ! do not activate it |
---|
[10025] | 215 | |
---|
| 216 | WRITE(*,*) '****** Bathymetry successfully created if no new topo ******' |
---|
| 217 | STOP |
---|
[2136] | 218 | ! |
---|
[10025] | 219 | ! ! ----------------------------------------------------- |
---|
| 220 | ELSE ! **** Second option : new topo file or partial steps |
---|
| 221 | ! ! ----------------------------------------------------- |
---|
[10381] | 222 | WRITE(*,*) '*** Second option : new topo or partial steps' |
---|
[10025] | 223 | |
---|
[2136] | 224 | ! read fine and coarse grids coordinates file |
---|
| 225 | status = Read_Coordinates(TRIM(parent_coordinate_file),G0) |
---|
[10381] | 226 | status = Read_Coordinates(TRIM(child_coordinates),G1,Pacifique) |
---|
[10025] | 227 | ! |
---|
| 228 | ! check error in size |
---|
| 229 | IF( imax > SIZE(G0%nav_lon,1) .OR. jmax > SIZE(G0%nav_lon,2) .OR. imax <= imin .OR. jmax <= jmin ) THEN |
---|
[2136] | 230 | WRITE(*,*) 'ERROR ***** bad child grid definition ...' |
---|
| 231 | WRITE(*,*) 'please check imin,jmin,imax,jmax,jpizoom,jpjzoom values' |
---|
| 232 | STOP |
---|
| 233 | ENDIF |
---|
| 234 | IF( SIZE(G1%nav_lon,1) .NE. nxfin .OR. SIZE(G1%nav_lon,2) .NE. nyfin ) THEN |
---|
| 235 | WRITE(*,*) 'ERROR ***** bad child coordinates file ...' |
---|
[10025] | 236 | WRITE(*,*) 'please check that child coordinates file has been created with the same namelist' |
---|
[2136] | 237 | STOP |
---|
| 238 | ENDIF |
---|
[10025] | 239 | ! |
---|
[10381] | 240 | ! === From here on: G0 is the grid associated with the new topography (as gebco or etopo) === |
---|
[10025] | 241 | ! |
---|
| 242 | ! read bathymetry data set => G0%bathy_meter |
---|
| 243 | IF( new_topo ) THEN ! read G0%bathy_meter (on a reduced grid) and G1 coordinates |
---|
[2136] | 244 | DEALLOCATE( G0%nav_lon, G0%nav_lat ) |
---|
[10381] | 245 | status = read_bathy_coord(TRIM(elevation_database),G0,G1,Pacifique) |
---|
[10025] | 246 | ELSE ! read G0%bathy_meter (on the global grid) |
---|
[10381] | 247 | IF( TRIM(parent_bathy_level) /= '' ) THEN |
---|
| 248 | status = Read_bathy_level(TRIM(parent_bathy_level),G0) |
---|
| 249 | CALL levels_to_meter(G0) |
---|
| 250 | ELSEIF( TRIM(parent_bathy_level) == '' .AND. TRIM(parent_bathy_meter) /= '') THEN |
---|
| 251 | status = read_bathy_meter(TRIM(parent_bathy_meter),G0) |
---|
| 252 | ENDIF |
---|
[2136] | 253 | IF(Pacifique) THEN |
---|
[10025] | 254 | WHERE(G0%nav_lon < 0.0001) G0%nav_lon = G0%nav_lon + 360. |
---|
[2136] | 255 | ENDIF |
---|
| 256 | ENDIF |
---|
| 257 | ! |
---|
[10025] | 258 | ! what type of interpolation for bathymetry |
---|
[2136] | 259 | IF( type_bathy_interp == 0 ) THEN |
---|
[10025] | 260 | WRITE(*,*) 'Interpolation of high resolution bathymetry on child grid: Arithmetic average ...' |
---|
[2136] | 261 | ELSE IF( type_bathy_interp == 1 ) THEN |
---|
[10025] | 262 | WRITE(*,*) 'Interpolation of high resolution bathymetry on child grid: Median average ...' |
---|
[2136] | 263 | ELSE IF( type_bathy_interp == 2 ) THEN |
---|
[10025] | 264 | WRITE(*,*) 'Interpolation of high resolution bathymetry on child grid: Bilinear interpolation ...' |
---|
[2136] | 265 | ELSE |
---|
[10025] | 266 | WRITE(*,*) 'bad value for type_bathy_interp variable ( must be 0, 1 or 2 )' |
---|
[2136] | 267 | STOP |
---|
| 268 | ENDIF |
---|
| 269 | ! |
---|
[10025] | 270 | ! 1st allocation of child grid bathy |
---|
| 271 | ALLOCATE(G1%bathy_meter(nxfin,nyfin)) |
---|
| 272 | G1%bathy_meter(:,:)=0. |
---|
[2136] | 273 | ! |
---|
[10025] | 274 | ! --------------------------------------------------------------------------------- |
---|
| 275 | ! === Bathymetry of the fine grid (step1) === |
---|
| 276 | ! --------------------------------------------------------------------------------- |
---|
| 277 | ! ==> It gives G1%bathy_meter from G0%bathy_meter |
---|
| 278 | ! --------------------------------------------------------------------------------- |
---|
| 279 | ! ! ----------------------------- |
---|
| 280 | IF( type_bathy_interp == 0 .OR. type_bathy_interp == 1 ) THEN ! arithmetic or median averages |
---|
| 281 | ! ! ----------------------------- |
---|
[2136] | 282 | ALLOCATE(trouble_points(nxfin,nyfin)) |
---|
[10025] | 283 | trouble_points(:,:) = 0 |
---|
[2136] | 284 | ! |
---|
[10027] | 285 | DO jj = 2, nyfin |
---|
| 286 | DO ji = 2, nxfin |
---|
[2136] | 287 | ! |
---|
[10025] | 288 | ! fine grid cell extension |
---|
[2136] | 289 | Cell_lonmin = MIN(G1%glamf(ji-1,jj-1),G1%glamf(ji,jj-1),G1%glamf(ji,jj),G1%glamf(ji-1,jj)) |
---|
| 290 | Cell_lonmax = MAX(G1%glamf(ji-1,jj-1),G1%glamf(ji,jj-1),G1%glamf(ji,jj),G1%glamf(ji-1,jj)) |
---|
| 291 | Cell_latmin = MIN(G1%gphif(ji-1,jj-1),G1%gphif(ji,jj-1),G1%gphif(ji,jj),G1%gphif(ji-1,jj)) |
---|
[9694] | 292 | Cell_latmax = MAX(G1%gphif(ji-1,jj-1),G1%gphif(ji,jj-1),G1%gphif(ji,jj),G1%gphif(ji-1,jj)) |
---|
[2136] | 293 | ! |
---|
[10025] | 294 | ! look for points in G0 (bathy dataset) contained in the fine grid cells |
---|
[2136] | 295 | iimin = 1 |
---|
| 296 | DO WHILE( G0%nav_lon(iimin,1) < Cell_lonmin ) |
---|
| 297 | iimin = iimin + 1 |
---|
| 298 | ENDDO |
---|
| 299 | ! |
---|
| 300 | jjmin = 1 |
---|
| 301 | DO WHILE( G0%nav_lat(iimin,jjmin) < Cell_latmin ) |
---|
| 302 | jjmin = jjmin + 1 |
---|
| 303 | ENDDO |
---|
| 304 | ! |
---|
| 305 | iimax = iimin |
---|
| 306 | DO WHILE( G0%nav_lon(iimax,1) <= Cell_lonmax ) |
---|
| 307 | iimax = iimax + 1 |
---|
[10025] | 308 | iimax = MIN( iimax,SIZE(G0%bathy_meter,1)) |
---|
[2136] | 309 | ENDDO |
---|
| 310 | ! |
---|
| 311 | jjmax = jjmin |
---|
| 312 | DO WHILE( G0%nav_lat(iimax,jjmax) <= Cell_latmax ) |
---|
| 313 | jjmax = jjmax + 1 |
---|
[10025] | 314 | jjmax = MIN( jjmax,SIZE(G0%bathy_meter,2)) |
---|
[2136] | 315 | ENDDO |
---|
| 316 | ! |
---|
[10025] | 317 | IF( ln_agrif_domain ) THEN |
---|
| 318 | iimax = iimax-1 |
---|
| 319 | jjmax = jjmax-1 |
---|
| 320 | ELSE |
---|
| 321 | iimax = MAX(iimin,iimax-1) |
---|
| 322 | jjmax = MAX(jjmin,jjmax-1) |
---|
| 323 | ENDIF |
---|
[2136] | 324 | ! |
---|
| 325 | iimin = MAX( iimin,1 ) |
---|
| 326 | jjmin = MAX( jjmin,1 ) |
---|
| 327 | iimax = MIN( iimax,SIZE(G0%bathy_meter,1)) |
---|
| 328 | jjmax = MIN( jjmax,SIZE(G0%bathy_meter,2)) |
---|
| 329 | |
---|
| 330 | nxhr = iimax - iimin + 1 |
---|
| 331 | nyhr = jjmax - jjmin + 1 |
---|
| 332 | |
---|
| 333 | IF( nxhr == 0 .OR. nyhr == 0 ) THEN |
---|
[10025] | 334 | ! |
---|
[2136] | 335 | trouble_points(ji,jj) = 1 |
---|
[10025] | 336 | ! |
---|
[2136] | 337 | ELSE |
---|
[10025] | 338 | ! |
---|
| 339 | ALLOCATE( vardep(nxhr,nyhr), mask_oce(nxhr,nyhr) ) |
---|
[2136] | 340 | vardep(:,:) = G0%bathy_meter(iimin:iimax,jjmin:jjmax) |
---|
[10025] | 341 | ! |
---|
| 342 | WHERE( vardep(:,:) .GT. 0. ) ; mask_oce = 1 ; |
---|
| 343 | ELSEWHERE ; mask_oce = 0 ; |
---|
| 344 | ENDWHERE |
---|
| 345 | ! |
---|
[10027] | 346 | nxyhr = nxhr*nyhr |
---|
| 347 | IF( SUM(mask_oce) < 0.5*(nxyhr) ) THEN ! if more than half of the points are on land then bathy fine = 0 |
---|
[2136] | 348 | G1%bathy_meter(ji,jj) = 0. |
---|
| 349 | ELSE |
---|
[10025] | 350 | IF( type_bathy_interp == 0 ) THEN ! Arithmetic average |
---|
| 351 | G1%bathy_meter(ji,jj) = SUM( vardep(:,:) * mask_oce(:,:) ) / SUM( mask_oce(:,:) ) |
---|
| 352 | ELSE ! Median average |
---|
[10027] | 353 | ALLOCATE(vardep1d(nxyhr)) |
---|
| 354 | vardep1d = RESHAPE(vardep,(/ nxyhr /) ) |
---|
| 355 | !!CALL ssort(vardep1d,nxyhr) |
---|
| 356 | CALL quicksort(vardep1d,1,nxyhr) |
---|
[2136] | 357 | ! |
---|
| 358 | ! Calculate median |
---|
[10027] | 359 | IF (MOD(nxyhr,2) .NE. 0) THEN |
---|
| 360 | G1%bathy_meter(ji,jj) = vardep1d( nxyhr/2 + 1 ) |
---|
[2136] | 361 | ELSE |
---|
[10027] | 362 | G1%bathy_meter(ji,jj) = 0.5 * ( vardep1d(nxyhr/2) + vardep1d(nxyhr/2+1) ) |
---|
[2136] | 363 | END IF |
---|
[9694] | 364 | DEALLOCATE(vardep1d) |
---|
[2136] | 365 | ENDIF |
---|
| 366 | ENDIF |
---|
| 367 | DEALLOCATE (mask_oce,vardep) |
---|
[10025] | 368 | ! |
---|
[2136] | 369 | ENDIF |
---|
| 370 | ENDDO |
---|
| 371 | ENDDO |
---|
| 372 | |
---|
| 373 | IF( SUM( trouble_points ) > 0 ) THEN |
---|
[10025] | 374 | PRINT*,'too much empty cells, proceed to bilinear interpolation' |
---|
[2136] | 375 | type_bathy_interp = 2 |
---|
| 376 | ENDIF |
---|
[10025] | 377 | |
---|
[2136] | 378 | ENDIF |
---|
[10025] | 379 | ! ! ----------------------------- |
---|
| 380 | IF( type_bathy_interp == 2) THEN ! Bilinear interpolation |
---|
| 381 | ! ! ----------------------------- |
---|
[2136] | 382 | identical_grids = .FALSE. |
---|
| 383 | |
---|
[10025] | 384 | IF( SIZE(G0%nav_lat,1) == SIZE(G1%nav_lat,1) .AND. SIZE(G0%nav_lat,2) == SIZE(G1%nav_lat,2) .AND. & |
---|
| 385 | SIZE(G0%nav_lon,1) == SIZE(G1%nav_lon,1) .AND. SIZE(G0%nav_lon,2) == SIZE(G1%nav_lon,2) ) THEN |
---|
[2136] | 386 | IF( MAXVAL( ABS(G0%nav_lat(:,:)- G1%nav_lat(:,:)) ) < 0.0001 .AND. & |
---|
[10025] | 387 | MAXVAL( ABS(G0%nav_lon(:,:)- G1%nav_lon(:,:)) ) < 0.0001 ) THEN |
---|
| 388 | PRINT*,'same grid between ',elevation_database,' and child domain' |
---|
[2136] | 389 | G1%bathy_meter = G0%bathy_meter |
---|
| 390 | identical_grids = .TRUE. |
---|
| 391 | ENDIF |
---|
| 392 | ENDIF |
---|
| 393 | |
---|
| 394 | IF( .NOT. identical_grids ) THEN |
---|
| 395 | |
---|
| 396 | ALLOCATE(masksrc(SIZE(G0%bathy_meter,1),SIZE(G0%bathy_meter,2))) |
---|
| 397 | ALLOCATE(bathy_test(nxfin,nyfin)) |
---|
| 398 | ! |
---|
[10025] | 399 | !Where(G0%bathy_meter.le.0.00001) |
---|
| 400 | ! masksrc = .false. |
---|
| 401 | !ElseWhere |
---|
| 402 | masksrc = .TRUE. |
---|
| 403 | !End where |
---|
[2136] | 404 | ! |
---|
| 405 | ! compute remapping matrix thanks to SCRIP package |
---|
[10025] | 406 | CALL get_remap_matrix(G0%nav_lat,G1%nav_lat,G0%nav_lon,G1%nav_lon,masksrc,matrix,src_add,dst_add) |
---|
| 407 | CALL make_remap(G0%bathy_meter,bathy_test,nxfin,nyfin,matrix,src_add,dst_add) |
---|
[2136] | 408 | ! |
---|
| 409 | G1%bathy_meter = bathy_test |
---|
| 410 | ! |
---|
| 411 | DEALLOCATE(masksrc) |
---|
| 412 | DEALLOCATE(bathy_test) |
---|
| 413 | |
---|
| 414 | ENDIF |
---|
| 415 | ! |
---|
| 416 | ENDIF |
---|
[10025] | 417 | ! --- |
---|
| 418 | ! At this stage bathymetry in meters has already been interpolated on fine grid |
---|
| 419 | ! => G1%bathy_meter(nxfin,nyfin) |
---|
[2136] | 420 | ! |
---|
[10025] | 421 | ! Also G0 was the grid from the new bathymetry data set (etopo, gebco...) and not the coarse grid |
---|
| 422 | ! --- |
---|
[2136] | 423 | ! |
---|
[10025] | 424 | ! --------------------------------------------------------------------------------- |
---|
| 425 | ! === Bathymetry of the fine grid (step2) === |
---|
| 426 | ! --------------------------------------------------------------------------------- |
---|
| 427 | ! ==> It gives an update of G1%bathy_meter and G1%bathy_level |
---|
| 428 | ! --------------------------------------------------------------------------------- |
---|
| 429 | ! From here on: G0 is the coarse grid |
---|
[2136] | 430 | ! |
---|
[10025] | 431 | ! Coarse grid bathymetry : G0%bathy_meter (on the global grid) |
---|
[10381] | 432 | IF( TRIM(parent_bathy_meter) /= '') THEN |
---|
| 433 | status = read_bathy_meter(TRIM(parent_bathy_meter),G0) |
---|
[10025] | 434 | ELSE |
---|
[10381] | 435 | status = Read_bathy_level(TRIM(parent_bathy_level),G0) |
---|
| 436 | CALL levels_to_meter(G0) |
---|
[10025] | 437 | ENDIF |
---|
| 438 | |
---|
| 439 | ! Coarse grid coordinatees : G0 coordinates |
---|
| 440 | DEALLOCATE(G0%nav_lat,G0%nav_lon) |
---|
| 441 | status = Read_coordinates(TRIM(parent_coordinate_file),G0) |
---|
| 442 | |
---|
| 443 | ! allocate temporary arrays |
---|
| 444 | IF (.NOT.ASSOCIATED(G0%gdepw_ps)) ALLOCATE(G0%gdepw_ps (SIZE(G0%bathy_meter,1),SIZE(G0%bathy_meter,2))) |
---|
| 445 | IF (.NOT.ASSOCIATED(G1%gdepw_ps)) ALLOCATE(G1%gdepw_ps (SIZE(G1%bathy_meter,1),SIZE(G1%bathy_meter,2))) |
---|
| 446 | IF (.NOT.ASSOCIATED(gdepw_ps_interp)) ALLOCATE(gdepw_ps_interp(SIZE(G1%bathy_meter,1),SIZE(G1%bathy_meter,2))) |
---|
| 447 | ! |
---|
| 448 | IF( ln_agrif_domain ) THEN |
---|
| 449 | boundary = npt_copy*irafx + nbghostcellsfine + 1 |
---|
| 450 | ELSE |
---|
| 451 | boundary = npt_copy*irafx |
---|
| 452 | ENDIF |
---|
[2136] | 453 | ! |
---|
[10025] | 454 | ! compute G0%gdepw_ps and G1%gdepw_ps |
---|
| 455 | CALL get_partial_steps(G0) |
---|
| 456 | CALL get_partial_steps(G1) |
---|
| 457 | CALL bathymetry_control(G0%Bathy_level) |
---|
| 458 | |
---|
| 459 | ! --------------------------------------- |
---|
| 460 | ! Bathymetry at the boundaries (npt_copy) |
---|
| 461 | ! --------------------------------------- |
---|
| 462 | ! 1st step: interpolate coarse bathy on the fine grid (using partial steps or not) |
---|
| 463 | IF( partial_steps .AND. ln_agrif_domain ) THEN |
---|
| 464 | CALL Check_interp(G0,gdepw_ps_interp) |
---|
| 465 | ELSE |
---|
| 466 | gdepw_ps_interp = 0. * G1%gdepw_ps |
---|
| 467 | !!CALL agrif_interp(G0%gdepw_ps,gdepw_ps_interp,'T') |
---|
| 468 | CALL init_constant_bathy(G0%gdepw_ps,gdepw_ps_interp) |
---|
| 469 | ENDIF |
---|
| 470 | |
---|
| 471 | IF (.NOT.ASSOCIATED(G1%wgt)) ALLOCATE(G1%wgt(SIZE(G1%bathy_meter,1),SIZE(G1%bathy_meter,2))) |
---|
| 472 | G1%wgt(:,:) = 0. |
---|
| 473 | IF ((.NOT.ASSOCIATED(G0%wgt)).AND.bathy_update) THEN |
---|
| 474 | ALLOCATE(G0%wgt(SIZE(G0%nav_lat,1),SIZE(G0%nav_lat,2))) |
---|
| 475 | G0%wgt(:,:) = 0. |
---|
| 476 | ENDIF |
---|
[2136] | 477 | ! |
---|
[10025] | 478 | ! 2nd step: copy parent bathymetry at the boundaries |
---|
| 479 | DO jj=1,nyfin ! West and East |
---|
| 480 | IF ( gdepw_ps_interp(nbghostcellsfine+1,jj) > 0. ) THEN |
---|
| 481 | G1%gdepw_ps(1:boundary,jj) = gdepw_ps_interp(1:boundary,jj) |
---|
| 482 | G1%wgt(1:boundary,jj) = 1. |
---|
| 483 | ELSE |
---|
| 484 | G1%gdepw_ps(1:nbghostcellsfine+1,jj)=0. |
---|
| 485 | ENDIF |
---|
[2136] | 486 | ! |
---|
[10025] | 487 | IF ( gdepw_ps_interp(nxfin-nbghostcellsfine,jj) > 0.) THEN |
---|
| 488 | G1%gdepw_ps(nxfin-boundary+1:nxfin,jj)=gdepw_ps_interp(nxfin-boundary+1:nxfin,jj) |
---|
| 489 | G1%wgt(nxfin-boundary+1:nxfin,jj) = 1. |
---|
| 490 | ELSE |
---|
| 491 | G1%gdepw_ps(nxfin-nbghostcellsfine:nxfin,jj) = 0. |
---|
| 492 | ENDIF |
---|
| 493 | END DO |
---|
| 494 | ! |
---|
| 495 | DO ji=1,nxfin ! South and North |
---|
| 496 | IF (gdepw_ps_interp(ji,nbghostcellsfine+1)>0.) THEN |
---|
| 497 | G1%gdepw_ps(ji,1:boundary) = gdepw_ps_interp(ji,1:boundary) |
---|
| 498 | G1%wgt(ji,1:boundary) = 1. |
---|
| 499 | ELSE |
---|
| 500 | G1%gdepw_ps(ji,1:nbghostcellsfine+1)=0. |
---|
| 501 | ENDIF |
---|
[2136] | 502 | ! |
---|
[10025] | 503 | IF (gdepw_ps_interp(ji,nyfin-nbghostcellsfine)>0.) THEN |
---|
| 504 | G1%gdepw_ps(ji,nyfin-boundary+1:nyfin)=gdepw_ps_interp(ji,nyfin-boundary+1:nyfin) |
---|
| 505 | G1%wgt(ji,nyfin-boundary+1:nyfin) = 1. |
---|
| 506 | ELSE |
---|
| 507 | G1%gdepw_ps(ji,nyfin-nbghostcellsfine:nyfin) = 0. |
---|
[9749] | 508 | ENDIF |
---|
[10025] | 509 | END DO |
---|
| 510 | ! |
---|
| 511 | !clem: recalculate interpolation everywhere before linear connection (useless to me) |
---|
| 512 | IF( partial_steps .AND. ln_agrif_domain ) THEN |
---|
| 513 | gdepw_ps_interp = 0. |
---|
| 514 | CALL Check_interp(G0,gdepw_ps_interp) |
---|
| 515 | ENDIF |
---|
| 516 | ! |
---|
| 517 | ! ------------------------------------------------------- |
---|
| 518 | ! Bathymetry between boundaries and interior (npt_connect) |
---|
| 519 | ! -------------------------------------------------------- |
---|
| 520 | ! Make linear connection (on npt_connect*irafx points) between the boundaries and the interior |
---|
| 521 | IF( ln_agrif_domain ) THEN |
---|
| 522 | boundary = (npt_copy + npt_connect)*irafx + nbghostcellsfine + 1 |
---|
| 523 | ELSE |
---|
| 524 | boundary = (npt_copy + npt_connect)*irafx |
---|
| 525 | ENDIF |
---|
[2136] | 526 | |
---|
[10025] | 527 | IF( npt_connect > 0 ) THEN |
---|
[9694] | 528 | WRITE(*,*) ' linear connection on ',npt_connect,'coarse grid points' |
---|
[2136] | 529 | |
---|
| 530 | wghts = 1. |
---|
[10025] | 531 | DO ji = boundary - npt_connect*irafx + 1 , boundary |
---|
| 532 | wghts = wghts - (1. / (npt_connect*irafx + 1. ) ) |
---|
| 533 | DO jj=1,nyfin |
---|
| 534 | IF (G1%gdepw_ps(nbghostcellsfine+1,jj) > 0.) G1%wgt(ji,jj) = MAX(wghts, G1%wgt(ji,jj)) |
---|
| 535 | END DO |
---|
| 536 | END DO |
---|
| 537 | |
---|
[9149] | 538 | wghts = 1. |
---|
[10025] | 539 | DO ji = nxfin - (boundary - npt_connect*irafx), nxfin - boundary +1 , -1 |
---|
| 540 | wghts = wghts - (1. / (npt_connect*irafx + 1. ) ) |
---|
| 541 | DO jj=1,nyfin |
---|
| 542 | IF (G1%gdepw_ps(nxfin-nbghostcellsfine,jj) > 0.) G1%wgt(ji,jj) = MAX(wghts, G1%wgt(ji,jj)) |
---|
| 543 | END DO |
---|
| 544 | END DO |
---|
[2136] | 545 | |
---|
| 546 | wghts = 1. |
---|
[10025] | 547 | DO jj = boundary - npt_connect*irafy + 1 , boundary |
---|
| 548 | wghts = wghts - (1. / (npt_connect*irafy + 1. ) ) |
---|
| 549 | DO ji=1,nxfin |
---|
| 550 | IF (G1%gdepw_ps(ji,nbghostcellsfine+1) > 0.) G1%wgt(ji,jj) = MAX(wghts, G1%wgt(ji,jj)) |
---|
| 551 | END DO |
---|
[9149] | 552 | END DO |
---|
| 553 | |
---|
[2136] | 554 | wghts = 1. |
---|
[10025] | 555 | DO jj = nyfin - (boundary - npt_connect*irafy) , nyfin - boundary +1, -1 |
---|
| 556 | wghts = wghts - (1. / (npt_connect*irafy + 1. ) ) |
---|
| 557 | DO ji=1,nxfin |
---|
| 558 | IF (G1%gdepw_ps(ji,nyfin-nbghostcellsfine) > 0.) G1%wgt(ji,jj) = MAX(wghts, G1%wgt(ji,jj)) |
---|
| 559 | END DO |
---|
[9149] | 560 | END DO |
---|
[9753] | 561 | IF (.NOT.identical_grids) THEN |
---|
| 562 | G1%gdepw_ps(:,:) = (1.-G1%wgt(:,:)) * G1%gdepw_ps(:,:) + gdepw_ps_interp(:,:)*G1%wgt(:,:) |
---|
| 563 | ENDIF |
---|
[2136] | 564 | |
---|
[10025] | 565 | ENDIF |
---|
[2136] | 566 | |
---|
[10025] | 567 | ! replace G1%bathy_meter by G1%gdepw_ps |
---|
| 568 | G1%bathy_meter = G1%gdepw_ps |
---|
| 569 | ! |
---|
| 570 | ! -------------------- |
---|
| 571 | ! Bathymetry smoothing |
---|
| 572 | ! -------------------- |
---|
| 573 | IF( smoothing .AND. (.NOT.identical_grids) ) THEN |
---|
| 574 | ! Chanut: smoothing everywhere then discard result in connection zone |
---|
| 575 | CALL smooth_topo(G1%gdepw_ps(:,:),nbiter) |
---|
| 576 | WHERE (G1%wgt(:,:)==0.) G1%bathy_meter(:,:) = G1%gdepw_ps(:,:) |
---|
| 577 | ELSE |
---|
| 578 | WRITE(*,*) 'No smoothing process only connection is carried out' |
---|
| 579 | ENDIF |
---|
| 580 | ! |
---|
| 581 | ! ------------------ |
---|
| 582 | ! Remove closed seas |
---|
| 583 | ! ------------------ |
---|
| 584 | IF (removeclosedseas) THEN |
---|
| 585 | ALLOCATE(bathy_test(nxfin,nyfin)) |
---|
| 586 | bathy_test=0. |
---|
| 587 | WHERE (G1%bathy_meter(1,:) .GT.0.) bathy_test(1,:)=1 |
---|
| 588 | WHERE (G1%bathy_meter(nxfin,:).GT.0.) bathy_test(nxfin,:)=1 |
---|
| 589 | WHERE (G1%bathy_meter(:,1) .GT.0.) bathy_test(:,1)=1 |
---|
| 590 | WHERE (G1%bathy_meter(:,nyfin).GT.0.) bathy_test(:,nyfin)=1 |
---|
[2136] | 591 | |
---|
[10025] | 592 | nbadd = 1 |
---|
| 593 | DO WHILE (nbadd.NE.0) |
---|
| 594 | nbadd = 0 |
---|
[10027] | 595 | DO jj=2,nyfin-1 |
---|
| 596 | DO ji=2,nxfin-1 |
---|
| 597 | IF (G1%bathy_meter(ji,jj).GT.0.) THEN |
---|
| 598 | IF (MAX(bathy_test(ji,jj+1),bathy_test(ji,jj-1),bathy_test(ji-1,jj),bathy_test(ji+1,jj)).EQ.1) THEN |
---|
| 599 | IF (bathy_test(ji,jj).NE.1.) nbadd = nbadd + 1 |
---|
| 600 | bathy_test(ji,jj)=1. |
---|
[10025] | 601 | ENDIF |
---|
[2136] | 602 | |
---|
[10025] | 603 | ENDIF |
---|
[2136] | 604 | ENDDO |
---|
| 605 | ENDDO |
---|
[10025] | 606 | ENDDO |
---|
| 607 | WHERE (bathy_test.EQ.0.) G1%bathy_meter = 0. |
---|
| 608 | DEALLOCATE(bathy_test) |
---|
| 609 | ENDIF |
---|
| 610 | ! |
---|
[10381] | 611 | IF( partial_steps ) THEN |
---|
| 612 | CALL get_partial_steps(G1) ! recompute bathy_level and gdepw_ps for G1 (and correct bathy_meter) |
---|
| 613 | ELSE |
---|
| 614 | CALL meter_to_levels(G1) ! convert bathymetry from meters to levels |
---|
[2136] | 615 | ENDIF |
---|
[10381] | 616 | ! |
---|
| 617 | ! update parent grid |
---|
| 618 | IF(bathy_update) CALL Update_Parent_Bathy( G0,G1 ) |
---|
| 619 | |
---|
| 620 | IF(bathy_update) status = Write_Bathy_meter(TRIM(parent_bathy_meter_updated),G0) |
---|
| 621 | IF(bathy_update) status = write_domcfg(TRIM(parent_domcfg_updated),G0) |
---|
| 622 | |
---|
| 623 | ! store interpolation result in output file |
---|
| 624 | IF( TRIM(parent_bathy_level) /= '' ) status = Write_Bathy_level(TRIM(child_level),G1) |
---|
| 625 | IF( TRIM(parent_bathy_meter) /= '' ) status = Write_Bathy_meter(TRIM(child_meter),G1) |
---|
| 626 | IF( TRIM(parent_domcfg_out) /= '' ) status = write_domcfg(TRIM(child_domcfg),G1) |
---|
| 627 | ! |
---|
| 628 | WRITE(*,*) '****** Bathymetry successfully created ******' |
---|
| 629 | STOP |
---|
[2136] | 630 | ENDIF |
---|
| 631 | ! |
---|
| 632 | END PROGRAM create_bathy |
---|
| 633 | |
---|
| 634 | |
---|