- Timestamp:
- 2019-12-16T18:42:48+01:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
utils/tools/NESTING/src/agrif_create_bathy.f90
r12253 r12259 84 84 STOP 85 85 ENDIF 86 ! 87 ! read fine and coarse grids coordinates file 88 status = Read_Coordinates(TRIM(parent_coordinate_file),G0) 89 status = Read_Coordinates(TRIM(child_coordinates),G1,Pacifique) 90 ! 91 ! check error in size 92 IF( imax > SIZE(G0%nav_lon,1) .OR. jmax > SIZE(G0%nav_lon,2) .OR. imax <= imin .OR. jmax <= jmin ) THEN 93 WRITE(*,*) 'ERROR ***** bad child grid definition ...' 94 WRITE(*,*) 'please check imin,jmin,imax,jmax,jpizoom,jpjzoom values' 95 STOP 96 ENDIF 97 IF( SIZE(G1%nav_lon,1) .NE. nxfin .OR. SIZE(G1%nav_lon,2) .NE. nyfin ) THEN 98 WRITE(*,*) 'ERROR ***** bad child coordinates file ...' 99 WRITE(*,*) 'please check that child coordinates file has been created with the same namelist' 100 STOP 101 ENDIF 102 ! 103 ! read bathymetry data set => G0%bathy_meter 104 IF( new_topo ) THEN ! read G0%bathy_meter (on a reduced grid) and G1 coordinates 105 DEALLOCATE( G0%nav_lon, G0%nav_lat ) 106 status = read_bathy_coord(TRIM(elevation_database),G0,G1,Pacifique) 107 CALL meter_to_levels(G0) 108 ELSE ! read G0%bathy_meter (on the global grid) 109 IF( TRIM(parent_bathy_level) /= '' ) THEN 110 status = Read_bathy_level(TRIM(parent_bathy_level),G0) 111 CALL levels_to_meter(G0) 112 ELSEIF( TRIM(parent_bathy_level) == '' .AND. TRIM(parent_bathy_meter) /= '') THEN 113 status = read_bathy_meter(TRIM(parent_bathy_meter),G0) 114 CALL meter_to_levels(G0) 115 ENDIF 116 ! change longitudes (from -180:180 to 0:360) 117 IF(Pacifique) THEN 118 WHERE(G0%nav_lon < 0.001) G0%nav_lon = G0%nav_lon + 360. 119 ENDIF 120 ENDIF 121 ! 122 ! 1st allocation of child grid bathy 123 ALLOCATE(G1%bathy_meter(nxfin,nyfin)) 124 G1%bathy_meter(:,:)=0. 125 126 ! check grids: if identical then do not interpolate 127 identical_grids = .FALSE. 128 129 IF( SIZE(G0%nav_lat,1) == SIZE(G1%nav_lat,1) .AND. SIZE(G0%nav_lat,2) == SIZE(G1%nav_lat,2) .AND. & 130 & SIZE(G0%nav_lon,1) == SIZE(G1%nav_lon,1) .AND. SIZE(G0%nav_lon,2) == SIZE(G1%nav_lon,2) ) THEN 131 IF( MAXVAL( ABS(G0%nav_lat(:,:)- G1%nav_lat(:,:)) ) < 0.0001 .AND. & 132 & MAXVAL( ABS(G0%nav_lon(:,:)- G1%nav_lon(:,:)) ) < 0.0001 ) THEN 133 WRITE(*,*) '' 134 WRITE(*,*) 'same grid between parent and child domains => NO INTERPOLATION' 135 WRITE(*,*) '' 136 G1%bathy_meter = G0%bathy_meter 137 identical_grids = .TRUE. 138 ENDIF 139 ENDIF 140 86 141 ! ! ------------------------------------------------------------------ 87 142 IF( .NOT.new_topo .AND. .NOT.partial_steps ) THEN ! **** First option : no new topo file & no partial steps … … 90 145 ! ! ------------------------------------------------------------------ 91 146 WRITE(*,*) '*** First option: no new topo file & no partial steps' 92 ! 93 ! read coarse grid coordinates 94 status = Read_Coordinates(TRIM(parent_coordinate_file),G0) 95 96 ! read coarse grid bathymetry 97 IF( TRIM(parent_bathy_level) /= '' ) THEN 98 status = Read_bathy_level(TRIM(parent_bathy_level),G0) 99 CALL levels_to_meter(G0) 100 ELSEIF( TRIM(parent_bathy_level) == '' .AND. TRIM(parent_bathy_meter) /= '') THEN 101 status = read_bathy_meter(TRIM(parent_bathy_meter),G0) 102 CALL meter_to_levels(G0) 103 ENDIF 104 ! 105 ! read fine grid coordinates 106 status = Read_Coordinates(TRIM(child_coordinates),G1,pacifique) 107 ! 108 ! stop if error in size 109 IF( imax > SIZE(G0%glamt,1) .OR. jmax > SIZE(G0%glamt,2) .OR. imax <= imin .OR. jmax <= jmin ) THEN 110 WRITE(*,*) 'ERROR ***** bad child grid definition ...' 111 WRITE(*,*) 'please check imin,jmin,imax,jmax,jpizoom,jpjzoom values' 112 STOP 113 ENDIF 114 IF( SIZE(G1%nav_lon,1) .NE. nxfin .OR. SIZE(G1%nav_lon,2) .NE. nyfin ) THEN 115 WRITE(*,*) 'ERROR ***** bad child coordinates file ...' 116 WRITE(*,*) 'please check that child coordinates file has been created with the same namelist' 117 STOP 118 ENDIF 119 ! 120 jpj = SIZE(G0%nav_lat,2) 121 jpi = SIZE(G0%nav_lat,1) 122 ! 123 ! create logical array masksrc 124 ALLOCATE(masksrc(jpi,jpj)) 125 WHERE(G0%bathy_level.LE.0) ; masksrc = .FALSE. ; 126 ELSEWHERE ; masksrc = .TRUE. ; 127 END WHERE 128 129 ! change longitudes (from -180:180 to 0:360) 130 IF ( Pacifique ) THEN 131 WHERE(G0%nav_lon < 0.001) G0%nav_lon = G0%nav_lon + 360. 132 ENDIF 133 ! 134 ALLOCATE(G1%bathy_meter(nxfin,nyfin)) 135 ! 136 ! compute remapping matrix thanks to SCRIP package 137 CALL get_remap_matrix(G0%nav_lat,G1%nav_lat,G0%nav_lon,G1%nav_lon,masksrc,matrix,src_add,dst_add) 138 CALL make_remap(G0%bathy_meter,G1%bathy_meter,nxfin,nyfin,matrix,src_add,dst_add) 139 DEALLOCATE(masksrc) 140 ! 147 ! 148 IF( .NOT. identical_grids ) THEN 149 ! 150 WRITE(*,*) 'Interpolation of high resolution bathymetry on child grid: Bilinear interpolation ...' 151 ! 152 ! create logical array masksrc 153 ALLOCATE(masksrc(SIZE(G0%bathy_meter,1),SIZE(G0%bathy_meter,2))) 154 WHERE(G0%bathy_level.LE.0) ; masksrc = .FALSE. ; 155 ELSEWHERE ; masksrc = .TRUE. ; 156 END WHERE 157 ! 158 ! compute remapping matrix thanks to SCRIP package 159 CALL get_remap_matrix(G0%nav_lat,G1%nav_lat,G0%nav_lon,G1%nav_lon,masksrc,matrix,src_add,dst_add) 160 CALL make_remap(G0%bathy_meter,G1%bathy_meter,nxfin,nyfin,matrix,src_add,dst_add) 161 DEALLOCATE(masksrc) 162 ! 163 ENDIF 164 141 165 ! compute constant bathymetry for Parent-Child bathymetry connection 142 166 CALL init_constant_bathy(G0%bathy_meter,bathy_fin_constant) … … 210 234 CALL levels_to_meter(G1) ! needed for domcfg 211 235 ! 236 ! update parent grid 237 IF(bathy_update) THEN 238 CALL Update_Parent_Bathy( G0,G1 ) 239 status = Write_Bathy_meter(TRIM(parent_bathy_meter_updated),G0) 240 status = write_domcfg(TRIM(parent_domcfg_updated),G0) 241 ENDIF 242 ! 212 243 ! store interpolation result in output file 213 status = Write_Bathy_level(TRIM(child_level),G1) 214 status = write_domcfg(TRIM(child_domcfg),G1) 244 IF( TRIM(parent_bathy_level) /= '' ) status = Write_Bathy_level(TRIM(child_level),G1) 245 IF( TRIM(parent_bathy_meter) /= '' ) status = Write_Bathy_meter(TRIM(child_meter),G1) 246 IF( TRIM(parent_domcfg_out) /= '' ) status = write_domcfg(TRIM(child_domcfg),G1) 215 247 !!status = write_domcfg(TRIM(parent_domcfg_out),G0) ! do not activate it 216 248 … … 222 254 ! ! ----------------------------------------------------- 223 255 WRITE(*,*) '*** Second option : new topo or partial steps' 224 225 ! read fine and coarse grids coordinates file 226 status = Read_Coordinates(TRIM(parent_coordinate_file),G0) 227 status = Read_Coordinates(TRIM(child_coordinates),G1,Pacifique) 228 ! 229 ! check error in size 230 IF( imax > SIZE(G0%nav_lon,1) .OR. jmax > SIZE(G0%nav_lon,2) .OR. imax <= imin .OR. jmax <= jmin ) THEN 231 WRITE(*,*) 'ERROR ***** bad child grid definition ...' 232 WRITE(*,*) 'please check imin,jmin,imax,jmax,jpizoom,jpjzoom values' 233 STOP 234 ENDIF 235 IF( SIZE(G1%nav_lon,1) .NE. nxfin .OR. SIZE(G1%nav_lon,2) .NE. nyfin ) THEN 236 WRITE(*,*) 'ERROR ***** bad child coordinates file ...' 237 WRITE(*,*) 'please check that child coordinates file has been created with the same namelist' 238 STOP 239 ENDIF 240 ! 241 ! === From here on: G0 is the grid associated with the new topography (as gebco or etopo) === 242 ! 243 ! read bathymetry data set => G0%bathy_meter 244 IF( new_topo ) THEN ! read G0%bathy_meter (on a reduced grid) and G1 coordinates 245 DEALLOCATE( G0%nav_lon, G0%nav_lat ) 246 status = read_bathy_coord(TRIM(elevation_database),G0,G1,Pacifique) 247 ELSE ! read G0%bathy_meter (on the global grid) 248 IF( TRIM(parent_bathy_level) /= '' ) THEN 249 status = Read_bathy_level(TRIM(parent_bathy_level),G0) 250 CALL levels_to_meter(G0) 251 ELSEIF( TRIM(parent_bathy_level) == '' .AND. TRIM(parent_bathy_meter) /= '') THEN 252 status = read_bathy_meter(TRIM(parent_bathy_meter),G0) 253 ENDIF 254 IF(Pacifique) THEN 255 WHERE(G0%nav_lon < 0.0001) G0%nav_lon = G0%nav_lon + 360. 256 ENDIF 257 ENDIF 258 ! 256 ! 257 ! === Here: G0 is the grid associated with the new topography (as gebco or etopo) === 258 ! 259 259 ! what type of interpolation for bathymetry 260 260 IF( type_bathy_interp == 0 ) THEN … … 269 269 ENDIF 270 270 ! 271 ! 1st allocation of child grid bathy272 ALLOCATE(G1%bathy_meter(nxfin,nyfin))273 G1%bathy_meter(:,:)=0.274 271 ! 275 272 ! --------------------------------------------------------------------------------- … … 278 275 ! ==> It gives G1%bathy_meter from G0%bathy_meter 279 276 ! --------------------------------------------------------------------------------- 280 identical_grids = .FALSE.281 282 IF( SIZE(G0%nav_lat,1) == SIZE(G1%nav_lat,1) .AND. SIZE(G0%nav_lat,2) == SIZE(G1%nav_lat,2) .AND. &283 & SIZE(G0%nav_lon,1) == SIZE(G1%nav_lon,1) .AND. SIZE(G0%nav_lon,2) == SIZE(G1%nav_lon,2) ) THEN284 IF( MAXVAL( ABS(G0%nav_lat(:,:)- G1%nav_lat(:,:)) ) < 0.0001 .AND. &285 & MAXVAL( ABS(G0%nav_lon(:,:)- G1%nav_lon(:,:)) ) < 0.0001 ) THEN286 WRITE(*,*) ''287 WRITE(*,*) 'same grid between parent and child domains => NO INTERPOLATION'288 WRITE(*,*) ''289 G1%bathy_meter = G0%bathy_meter290 identical_grids = .TRUE.291 ENDIF292 ENDIF293 277 294 278 IF( .NOT. identical_grids ) THEN … … 392 376 ENDIF 393 377 378 DEALLOCATE(trouble_points) 379 394 380 ENDIF 395 381 ! ! ----------------------------- … … 400 386 ALLOCATE(bathy_test(nxfin,nyfin)) 401 387 ! 402 !Where(G0%bathy_meter.le.0.00001) 403 ! masksrc = .false. 404 !ElseWhere 405 masksrc = .TRUE. 406 !End where 388 WHERE(G0%bathy_level.LE.0) ; masksrc = .FALSE. ; 389 ELSEWHERE ; masksrc = .TRUE. ; 390 END WHERE 407 391 ! 408 392 ! compute remapping matrix thanks to SCRIP package … … 446 430 ! allocate temporary arrays 447 431 IF (.NOT.ASSOCIATED(G0%gdepw_ps)) ALLOCATE(G0%gdepw_ps (SIZE(G0%bathy_meter,1),SIZE(G0%bathy_meter,2))) 448 IF (.NOT.ASSOCIATED(G1%gdepw_ps)) ALLOCATE(G1%gdepw_ps (SIZE(G1%bathy_meter,1),SIZE(G1%bathy_meter,2))) 432 IF (.NOT.ASSOCIATED(G1%gdepw_ps)) ALLOCATE(G1%gdepw_ps (SIZE(G1%bathy_meter,1),SIZE(G1%bathy_meter,2))) 449 433 IF (.NOT.ASSOCIATED(gdepw_ps_interp)) ALLOCATE(gdepw_ps_interp(SIZE(G1%bathy_meter,1),SIZE(G1%bathy_meter,2))) 450 434 ! … … 619 603 ! 620 604 ! update parent grid 621 IF(bathy_update) CALL Update_Parent_Bathy( G0,G1 ) 622 623 IF(bathy_update) status = Write_Bathy_meter(TRIM(parent_bathy_meter_updated),G0) 624 IF(bathy_update) status = write_domcfg(TRIM(parent_domcfg_updated),G0) 625 605 IF(bathy_update) THEN 606 CALL Update_Parent_Bathy( G0,G1 ) 607 status = Write_Bathy_meter(TRIM(parent_bathy_meter_updated),G0) 608 status = write_domcfg(TRIM(parent_domcfg_updated),G0) 609 ENDIF 610 ! 626 611 ! store interpolation result in output file 627 612 IF( TRIM(parent_bathy_level) /= '' ) status = Write_Bathy_level(TRIM(child_level),G1)
Note: See TracChangeset
for help on using the changeset viewer.