Changeset 12412
- Timestamp:
- 2020-02-19T18:48:10+01:00 (4 years ago)
- Location:
- utils/tools_MERGE_2019
- Files:
-
- 8 edited
-
NESTING/src/agrif_connect_topo.f90 (modified) (1 diff)
-
NESTING/src/agrif_create_bathy.f90 (modified) (2 diffs)
-
NESTING/src/agrif_interpolation.f90 (modified) (2 diffs)
-
NESTING/src/agrif_partial_steps.f90 (modified) (5 diffs)
-
NESTING/src/agrif_readwrite.f90 (modified) (4 diffs)
-
NESTING/src/agrif_types.f90 (modified) (1 diff)
-
SIREN/src/create_boundary.F90 (modified) (1 diff)
-
SIREN/src/docsrc/5_changeLog.md (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
utils/tools_MERGE_2019/NESTING/src/agrif_connect_topo.f90
r10025 r12412 251 251 ! 252 252 ! 253 DO j = 1,jpj254 DO i = 1,jpi253 DO jj = 1,jpj 254 DO ji = 1,jpi 255 255 ! 256 IF (Grid%bathy_meter( i,j) .EQ. 0.0 ) THEN257 Grid%bathy_level( i,j)=0.256 IF (Grid%bathy_meter(ji,jj) .EQ. 0.0 ) THEN 257 Grid%bathy_level(ji,jj)=0 258 258 ELSE 259 259 ! 260 k1= 4260 k1=2 ! clem: minimum levels = 2 ??? 261 261 DO WHILE (k1 .LT. (N-1)) 262 IF ((Grid%bathy_meter( i,j).GE.gdepw(k1)) &263 .AND.(Grid%bathy_meter( i,j).LE.gdepw(k1+1))) EXIT262 IF ((Grid%bathy_meter(ji,jj).GE.gdepw(k1)) & 263 .AND.(Grid%bathy_meter(ji,jj).LE.gdepw(k1+1))) EXIT 264 264 k1=k1+1 265 265 END DO 266 Grid%bathy_level( i,j)=k1266 Grid%bathy_level(ji,jj)=k1 267 267 ! 268 268 ENDIF -
utils/tools_MERGE_2019/NESTING/src/agrif_create_bathy.f90
r10381 r12412 84 84 STOP 85 85 ENDIF 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' 91 ! 92 ! read coarse grid coordinates 93 status = Read_Coordinates(TRIM(parent_coordinate_file),G0) 94 95 ! read coarse grid bathymetry 96 IF( TRIM(parent_bathy_level) /= '' ) THEN 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 ELSE ! read G0%bathy_meter (on the global grid) 108 IF( TRIM(parent_bathy_meter) /= '') THEN 109 status = read_bathy_meter(TRIM(parent_bathy_meter),G0) 110 ELSE 97 111 status = Read_bathy_level(TRIM(parent_bathy_level),G0) 98 112 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 103 ! 104 ! read fine grid coordinates 105 status = Read_Coordinates(TRIM(child_coordinates),G1,pacifique) 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 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 ...' 115 WRITE(*,*) 'please check that child coordinates file has been created with the same namelist' 116 STOP 117 ENDIF 118 ! 119 jpj = SIZE(G0%nav_lat,2) 120 jpi = SIZE(G0%nav_lat,1) 121 ! 122 ! create logical array masksrc 123 ALLOCATE(masksrc(jpi,jpj)) 124 WHERE(G0%bathy_level.LE.0) ; masksrc = .FALSE. ; 125 ELSEWHERE ; masksrc = .TRUE. ; 126 END WHERE 127 113 ENDIF 128 114 ! change longitudes (from -180:180 to 0:360) 129 IF ( Pacifique) THEN115 IF(Pacifique) THEN 130 116 WHERE(G0%nav_lon < 0.001) G0%nav_lon = G0%nav_lon + 360. 131 117 ENDIF 132 ! 133 ALLOCATE(G1%bathy_meter(nxfin,nyfin)) 134 ! 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) 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 ! 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 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) 153 DEALLOCATE(bathy_fin_constant) 154 ! 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 165 CALL meter_to_levels(G1) 166 G1%bathy_level=NINT(G1%bathy_level) 167 ! 168 ! Check errors in bathy 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) 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 178 ! 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 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. 198 ENDIF 199 ENDIF 200 ENDDO 201 ENDDO 202 ENDDO 203 204 WHERE (bathy_test.EQ.0.) G1%bathy_level = 0. 205 206 DEALLOCATE(bathy_test) 207 ENDIF 208 ! 209 CALL levels_to_meter(G1) ! needed for domcfg 210 ! 211 ! store interpolation result in output file 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 215 216 WRITE(*,*) '****** Bathymetry successfully created if no new topo ******' 217 STOP 218 ! 219 ! ! ----------------------------------------------------- 220 ELSE ! **** Second option : new topo file or partial steps 221 ! ! ----------------------------------------------------- 222 WRITE(*,*) '*** Second option : new topo or partial steps' 223 224 ! read fine and coarse grids coordinates file 225 status = Read_Coordinates(TRIM(parent_coordinate_file),G0) 226 status = Read_Coordinates(TRIM(child_coordinates),G1,Pacifique) 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 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 ...' 236 WRITE(*,*) 'please check that child coordinates file has been created with the same namelist' 237 STOP 238 ENDIF 239 ! 240 ! === From here on: G0 is the grid associated with the new topography (as gebco or etopo) === 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 244 DEALLOCATE( G0%nav_lon, G0%nav_lat ) 245 status = read_bathy_coord(TRIM(elevation_database),G0,G1,Pacifique) 246 ELSE ! read G0%bathy_meter (on the global grid) 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 253 IF(Pacifique) THEN 254 WHERE(G0%nav_lon < 0.0001) G0%nav_lon = G0%nav_lon + 360. 255 ENDIF 256 ENDIF 257 ! 258 ! what type of interpolation for bathymetry 259 IF( type_bathy_interp == 0 ) THEN 260 WRITE(*,*) 'Interpolation of high resolution bathymetry on child grid: Arithmetic average ...' 261 ELSE IF( type_bathy_interp == 1 ) THEN 262 WRITE(*,*) 'Interpolation of high resolution bathymetry on child grid: Median average ...' 263 ELSE IF( type_bathy_interp == 2 ) THEN 264 WRITE(*,*) 'Interpolation of high resolution bathymetry on child grid: Bilinear interpolation ...' 265 ELSE 266 WRITE(*,*) 'bad value for type_bathy_interp variable ( must be 0, 1 or 2 )' 267 STOP 268 ENDIF 269 ! 270 ! 1st allocation of child grid bathy 271 ALLOCATE(G1%bathy_meter(nxfin,nyfin)) 272 G1%bathy_meter(:,:)=0. 273 ! 274 ! --------------------------------------------------------------------------------- 275 ! === Bathymetry of the fine grid (step1) === 276 ! --------------------------------------------------------------------------------- 277 ! ==> It gives G1%bathy_meter from G0%bathy_meter 278 ! --------------------------------------------------------------------------------- 118 ENDIF 119 ! 120 ! 1st allocation of child grid bathy 121 ALLOCATE(G1%bathy_meter(nxfin,nyfin)) 122 G1%bathy_meter(:,:)=0. 123 124 ! check grids: if identical then do not interpolate 125 identical_grids = .FALSE. 126 127 IF( SIZE(G0%nav_lat,1) == SIZE(G1%nav_lat,1) .AND. SIZE(G0%nav_lat,2) == SIZE(G1%nav_lat,2) .AND. & 128 & SIZE(G0%nav_lon,1) == SIZE(G1%nav_lon,1) .AND. SIZE(G0%nav_lon,2) == SIZE(G1%nav_lon,2) ) THEN 129 IF( MAXVAL( ABS(G0%nav_lat(:,:)- G1%nav_lat(:,:)) ) < 0.0001 .AND. & 130 & MAXVAL( ABS(G0%nav_lon(:,:)- G1%nav_lon(:,:)) ) < 0.0001 ) THEN 131 WRITE(*,*) '' 132 WRITE(*,*) 'same grid between parent and child domains => NO INTERPOLATION' 133 WRITE(*,*) '' 134 G1%bathy_meter = G0%bathy_meter 135 identical_grids = .TRUE. 136 ENDIF 137 ENDIF 138 139 IF( .NOT.new_topo ) type_bathy_interp = 2 ! only one which works 140 ! 141 ! 142 ! what type of interpolation for bathymetry 143 IF( type_bathy_interp == 0 ) THEN 144 WRITE(*,*) 'Interpolation of high resolution bathymetry on child grid: Arithmetic average ...' 145 ELSE IF( type_bathy_interp == 1 ) THEN 146 WRITE(*,*) 'Interpolation of high resolution bathymetry on child grid: Median average ...' 147 ELSE IF( type_bathy_interp == 2 ) THEN 148 WRITE(*,*) 'Interpolation of high resolution bathymetry on child grid: Bilinear interpolation ...' 149 ELSE 150 WRITE(*,*) 'bad value for type_bathy_interp variable ( must be 0, 1 or 2 )' 151 STOP 152 ENDIF 153 ! 154 ! 155 ! --------------------------------------------------------------------------------- 156 ! === Bathymetry of the fine grid (step1) === 157 ! --------------------------------------------------------------------------------- 158 ! ==> It gives G1%bathy_meter from G0%bathy_meter 159 ! --------------------------------------------------------------------------------- 160 161 ! === Here: G0 is the grid associated with the new topography (as gebco or etopo) === 162 163 IF( .NOT. identical_grids ) THEN 279 164 ! ! ----------------------------- 280 165 IF( type_bathy_interp == 0 .OR. type_bathy_interp == 1 ) THEN ! arithmetic or median averages … … 375 260 type_bathy_interp = 2 376 261 ENDIF 377 262 263 DEALLOCATE(trouble_points) 264 378 265 ENDIF 379 266 ! ! ----------------------------- 380 267 IF( type_bathy_interp == 2) THEN ! Bilinear interpolation 381 268 ! ! ----------------------------- 382 identical_grids = .FALSE. 383 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 386 IF( MAXVAL( ABS(G0%nav_lat(:,:)- G1%nav_lat(:,:)) ) < 0.0001 .AND. & 387 MAXVAL( ABS(G0%nav_lon(:,:)- G1%nav_lon(:,:)) ) < 0.0001 ) THEN 388 PRINT*,'same grid between ',elevation_database,' and child domain' 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 ! 399 !Where(G0%bathy_meter.le.0.00001) 400 ! masksrc = .false. 401 !ElseWhere 402 masksrc = .TRUE. 403 !End where 404 ! 405 ! compute remapping matrix thanks to SCRIP package 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) 408 ! 409 G1%bathy_meter = bathy_test 410 ! 411 DEALLOCATE(masksrc) 412 DEALLOCATE(bathy_test) 413 414 ENDIF 269 270 ALLOCATE(masksrc(SIZE(G0%bathy_meter,1),SIZE(G0%bathy_meter,2))) 271 ALLOCATE(bathy_test(nxfin,nyfin)) 272 ! 273 WHERE(G0%bathy_meter.LE.0) ; masksrc = .FALSE. ; 274 ELSEWHERE ; masksrc = .TRUE. ; 275 END WHERE 415 276 ! 416 ENDIF 417 ! --- 418 ! At this stage bathymetry in meters has already been interpolated on fine grid 419 ! => G1%bathy_meter(nxfin,nyfin) 277 ! compute remapping matrix thanks to SCRIP package 278 CALL get_remap_matrix(G0%nav_lat,G1%nav_lat,G0%nav_lon,G1%nav_lon,masksrc,matrix,src_add,dst_add) 279 CALL make_remap(G0%bathy_meter,bathy_test,nxfin,nyfin,matrix,src_add,dst_add) 280 ! 281 G1%bathy_meter = bathy_test 282 ! 283 DEALLOCATE(masksrc) 284 DEALLOCATE(bathy_test) 285 286 ENDIF 287 ! 288 ENDIF ! not identical grids 289 ! --- 290 ! At this stage bathymetry in meters has already been interpolated on fine grid 291 ! => G1%bathy_meter(nxfin,nyfin) 292 ! 293 ! Also G0 was the grid from the new bathymetry data set (etopo, gebco...) and not the coarse grid 294 ! --- 295 ! 296 ! --------------------------------------------------------------------------------- 297 ! === Bathymetry of the fine grid (step2) === 298 ! --------------------------------------------------------------------------------- 299 ! ==> It gives an update of G1%bathy_meter and G1%bathy_level 300 ! --------------------------------------------------------------------------------- 301 ! From here on: G0 is the coarse grid 302 ! 303 ! Coarse grid bathymetry : G0%bathy_meter (on the global grid) 304 IF( TRIM(parent_bathy_meter) /= '') THEN 305 status = read_bathy_meter(TRIM(parent_bathy_meter),G0) 306 ELSE 307 status = Read_bathy_level(TRIM(parent_bathy_level),G0) 308 CALL levels_to_meter(G0) 309 ENDIF 310 311 ! Coarse grid coordinatees : G0 coordinates 312 DEALLOCATE(G0%nav_lat,G0%nav_lon) 313 status = Read_coordinates(TRIM(parent_coordinate_file),G0) 314 315 ! allocate temporary arrays 316 IF (.NOT.ASSOCIATED(G0%gdepw_ps)) ALLOCATE(G0%gdepw_ps (SIZE(G0%bathy_meter,1),SIZE(G0%bathy_meter,2))) 317 IF (.NOT.ASSOCIATED(G1%gdepw_ps)) ALLOCATE(G1%gdepw_ps (SIZE(G1%bathy_meter,1),SIZE(G1%bathy_meter,2))) 318 IF (.NOT.ASSOCIATED(gdepw_ps_interp)) ALLOCATE(gdepw_ps_interp(SIZE(G1%bathy_meter,1),SIZE(G1%bathy_meter,2))) 319 ! 320 IF( ln_agrif_domain ) THEN 321 boundary = npt_copy*irafx + nbghostcellsfine + 1 322 ELSE 323 boundary = npt_copy*irafx 324 ENDIF 325 ! 326 ! compute G0%gdepw_ps and G1%gdepw_ps 327 CALL get_partial_steps(G0) 328 CALL get_partial_steps(G1) 329 CALL bathymetry_control(G0%Bathy_level) 330 331 ! --------------------------------------- 332 ! Bathymetry at the boundaries (npt_copy) 333 ! --------------------------------------- 334 ! 1st step: interpolate coarse bathy on the fine grid (using partial steps or not) 335 IF( ln_agrif_domain ) THEN 336 CALL Check_interp(G0,gdepw_ps_interp) 337 ELSE 338 gdepw_ps_interp = 0. * G1%gdepw_ps 339 !!CALL agrif_interp(G0%gdepw_ps,gdepw_ps_interp,'T') 340 CALL init_constant_bathy(G0%gdepw_ps,gdepw_ps_interp) 341 ENDIF 342 343 IF (.NOT.ASSOCIATED(G1%wgt)) ALLOCATE(G1%wgt(SIZE(G1%bathy_meter,1),SIZE(G1%bathy_meter,2))) 344 G1%wgt(:,:) = 0. 345 IF ((.NOT.ASSOCIATED(G0%wgt)).AND.bathy_update) THEN 346 ALLOCATE(G0%wgt(SIZE(G0%nav_lat,1),SIZE(G0%nav_lat,2))) 347 G0%wgt(:,:) = 0. 348 ENDIF 349 ! 350 !!$ IF( new_topo ) THEN ! clem: no, do it even when there is no new topo 351 ! 2nd step: copy parent bathymetry at the boundaries 352 DO jj=1,nyfin ! West and East 353 IF ( gdepw_ps_interp(nbghostcellsfine+1,jj) > 0. ) THEN 354 G1%gdepw_ps(1:boundary,jj) = gdepw_ps_interp(1:boundary,jj) 355 G1%wgt(1:boundary,jj) = 1. 356 ELSE 357 G1%gdepw_ps(1:nbghostcellsfine+1,jj)=0. 358 ENDIF 420 359 ! 421 ! Also G0 was the grid from the new bathymetry data set (etopo, gebco...) and not the coarse grid 422 ! --- 360 IF ( gdepw_ps_interp(nxfin-nbghostcellsfine,jj) > 0.) THEN 361 G1%gdepw_ps(nxfin-boundary+1:nxfin,jj)=gdepw_ps_interp(nxfin-boundary+1:nxfin,jj) 362 G1%wgt(nxfin-boundary+1:nxfin,jj) = 1. 363 ELSE 364 G1%gdepw_ps(nxfin-nbghostcellsfine:nxfin,jj) = 0. 365 ENDIF 366 END DO 367 ! 368 DO ji=1,nxfin ! South and North 369 IF (gdepw_ps_interp(ji,nbghostcellsfine+1)>0.) THEN 370 G1%gdepw_ps(ji,1:boundary) = gdepw_ps_interp(ji,1:boundary) 371 G1%wgt(ji,1:boundary) = 1. 372 ELSE 373 G1%gdepw_ps(ji,1:nbghostcellsfine+1)=0. 374 ENDIF 423 375 ! 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 430 ! 431 ! Coarse grid bathymetry : G0%bathy_meter (on the global grid) 432 IF( TRIM(parent_bathy_meter) /= '') THEN 433 status = read_bathy_meter(TRIM(parent_bathy_meter),G0) 376 IF (gdepw_ps_interp(ji,nyfin-nbghostcellsfine)>0.) THEN 377 G1%gdepw_ps(ji,nyfin-boundary+1:nyfin)=gdepw_ps_interp(ji,nyfin-boundary+1:nyfin) 378 G1%wgt(ji,nyfin-boundary+1:nyfin) = 1. 434 379 ELSE 435 status = Read_bathy_level(TRIM(parent_bathy_level),G0) 436 CALL levels_to_meter(G0) 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 453 ! 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 477 ! 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 486 ! 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 380 G1%gdepw_ps(ji,nyfin-nbghostcellsfine:nyfin) = 0. 381 ENDIF 382 END DO 383 ! 384 !clem: recalculate interpolation everywhere before linear connection (useless to me??) 385 IF( ln_agrif_domain ) THEN 386 gdepw_ps_interp = 0. 387 CALL Check_interp(G0,gdepw_ps_interp) 388 ENDIF 389 ! 390 ! ------------------------------------------------------- 391 ! Bathymetry between boundaries and interior (npt_connect) 392 ! -------------------------------------------------------- 393 ! Make linear connection (on npt_connect*irafx points) between the boundaries and the interior 394 IF( ln_agrif_domain ) THEN 395 boundary = (npt_copy + npt_connect)*irafx + nbghostcellsfine + 1 396 ELSE 397 boundary = (npt_copy + npt_connect)*irafx 398 ENDIF 399 400 IF( npt_connect > 0 ) THEN 401 WRITE(*,*) ' linear connection on ',npt_connect,'coarse grid points' 402 403 wghts = 1. 404 DO ji = boundary - npt_connect*irafx + 1 , boundary 405 wghts = wghts - (1. / (npt_connect*irafx + 1. ) ) 406 DO jj=1,nyfin 407 IF (G1%gdepw_ps(nbghostcellsfine+1,jj) > 0.) G1%wgt(ji,jj) = MAX(wghts, G1%wgt(ji,jj)) 408 END DO 493 409 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 502 ! 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. 508 ENDIF 410 411 wghts = 1. 412 DO ji = nxfin - (boundary - npt_connect*irafx), nxfin - boundary +1 , -1 413 wghts = wghts - (1. / (npt_connect*irafx + 1. ) ) 414 DO jj=1,nyfin 415 IF (G1%gdepw_ps(nxfin-nbghostcellsfine,jj) > 0.) G1%wgt(ji,jj) = MAX(wghts, G1%wgt(ji,jj)) 416 END DO 509 417 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 526 527 IF( npt_connect > 0 ) THEN 528 WRITE(*,*) ' linear connection on ',npt_connect,'coarse grid points' 529 530 wghts = 1. 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 418 419 wghts = 1. 420 DO jj = boundary - npt_connect*irafy + 1 , boundary 421 wghts = wghts - (1. / (npt_connect*irafy + 1. ) ) 422 DO ji=1,nxfin 423 IF (G1%gdepw_ps(ji,nbghostcellsfine+1) > 0.) G1%wgt(ji,jj) = MAX(wghts, G1%wgt(ji,jj)) 536 424 END DO 537 538 wghts = 1. 539 DO ji = nxfin - (boundary - npt_connect*irafx), nxfin - boundary +1 , -1540 wghts = wghts - (1. / (npt_connect*irafx + 1. ) )541 DO jj=1,nyfin542 IF (G1%gdepw_ps(nxfin-nbghostcellsfine,jj) > 0.) G1%wgt(ji,jj) = MAX(wghts, G1%wgt(ji,jj))543 END DO425 END DO 426 427 wghts = 1. 428 DO jj = nyfin - (boundary - npt_connect*irafy) , nyfin - boundary +1, -1 429 wghts = wghts - (1. / (npt_connect*irafy + 1. ) ) 430 DO ji=1,nxfin 431 IF (G1%gdepw_ps(ji,nyfin-nbghostcellsfine) > 0.) G1%wgt(ji,jj) = MAX(wghts, G1%wgt(ji,jj)) 544 432 END DO 545 546 wghts = 1. 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 552 END DO 553 554 wghts = 1. 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 560 END DO 561 IF (.NOT.identical_grids) THEN 562 G1%gdepw_ps(:,:) = (1.-G1%wgt(:,:)) * G1%gdepw_ps(:,:) + gdepw_ps_interp(:,:)*G1%wgt(:,:) 563 ENDIF 564 565 ENDIF 566 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 591 592 nbadd = 1 593 DO WHILE (nbadd.NE.0) 594 nbadd = 0 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. 601 ENDIF 602 433 END DO 434 IF (.NOT.identical_grids) THEN 435 G1%gdepw_ps(:,:) = (1.-G1%wgt(:,:)) * G1%gdepw_ps(:,:) + gdepw_ps_interp(:,:)*G1%wgt(:,:) 436 ENDIF 437 438 ENDIF 439 !!$ ENDIF 440 441 ! replace G1%bathy_meter by G1%gdepw_ps 442 G1%bathy_meter = G1%gdepw_ps 443 ! 444 ! -------------------- 445 ! Bathymetry smoothing 446 ! -------------------- 447 IF( smoothing .AND. (.NOT.identical_grids) ) THEN 448 ! Chanut: smoothing everywhere then discard result in connection zone 449 CALL smooth_topo(G1%gdepw_ps(:,:),nbiter) 450 WHERE (G1%wgt(:,:)==0.) G1%bathy_meter(:,:) = G1%gdepw_ps(:,:) 451 ELSE 452 WRITE(*,*) 'No smoothing process only connection is carried out' 453 ENDIF 454 ! 455 ! ------------------ 456 ! Remove closed seas 457 ! ------------------ 458 IF (removeclosedseas) THEN 459 ALLOCATE(bathy_test(nxfin,nyfin)) 460 bathy_test=0. 461 WHERE (G1%bathy_meter(1,:) .GT.0.) bathy_test(1,:)=1 462 WHERE (G1%bathy_meter(nxfin,:).GT.0.) bathy_test(nxfin,:)=1 463 WHERE (G1%bathy_meter(:,1) .GT.0.) bathy_test(:,1)=1 464 WHERE (G1%bathy_meter(:,nyfin).GT.0.) bathy_test(:,nyfin)=1 465 466 nbadd = 1 467 DO WHILE (nbadd.NE.0) 468 nbadd = 0 469 DO jj=2,nyfin-1 470 DO ji=2,nxfin-1 471 IF (G1%bathy_meter(ji,jj).GT.0.) THEN 472 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 473 IF (bathy_test(ji,jj).NE.1.) nbadd = nbadd + 1 474 bathy_test(ji,jj)=1. 603 475 ENDIF 604 ENDDO 476 477 ENDIF 605 478 ENDDO 606 479 ENDDO 607 WHERE (bathy_test.EQ.0.) G1%bathy_meter = 0. 608 DEALLOCATE(bathy_test) 609 ENDIF 610 ! 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 615 ENDIF 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 630 ENDIF 480 ENDDO 481 WHERE (bathy_test.EQ.0.) G1%bathy_meter = 0. 482 DEALLOCATE(bathy_test) 483 ENDIF 484 ! 485 CALL get_partial_steps(G1) ! recompute bathy_level and gdepw_ps for G1 (and correct bathy_meter) 486 ! 487 ! update parent grid 488 IF(bathy_update) THEN 489 CALL Update_Parent_Bathy( G0,G1 ) 490 status = Write_Bathy_meter(TRIM(parent_bathy_meter_updated),G0) 491 status = write_domcfg(TRIM(parent_domcfg_updated),G0) 492 ENDIF 493 ! 494 ! store interpolation result in output file 495 IF( TRIM(parent_bathy_level) /= '' ) status = Write_Bathy_level(TRIM(child_level),G1) 496 IF( TRIM(parent_bathy_meter) /= '' ) status = Write_Bathy_meter(TRIM(child_meter),G1) 497 IF( TRIM(parent_domcfg_out) /= '' ) status = write_domcfg(TRIM(child_domcfg),G1) 498 ! 499 WRITE(*,*) '****** Bathymetry successfully created ******' 500 STOP 631 501 ! 632 502 END PROGRAM create_bathy -
utils/tools_MERGE_2019/NESTING/src/agrif_interpolation.f90
r10381 r12412 198 198 CHARACTER(*) :: typevar 199 199 200 INTEGER :: nxf,nyf, zx,zy201 INTEGER :: ji,jj,jif,jjf,jic,jjc,j decx,jdecy200 INTEGER :: nxf,nyf,nxc,nyc,zx,zy 201 INTEGER :: ji,jj,jif,jjf,jic,jjc,jic1,jjc1,jdecx,jdecy 202 202 REAL*8 :: Ax, Bx, Ay, By 203 203 204 204 nxf = SIZE(tabout,1) 205 205 nyf = SIZE(tabout,2) 206 207 nxc = SIZE(tabin,1) 208 nyc = SIZE(tabin,2) 206 209 207 210 SELECT CASE(typevar) … … 248 251 Ay = 1. - By 249 252 250 tabout(ji,jj) = ( Bx * tabin(jic+1,jjc ) + Ax * tabin(jic,jjc ) ) * Ay + & 251 & ( Bx * tabin(jic+1,jjc+1) + Ax * tabin(jic,jjc+1) ) * By 253 jic1 = MIN( nxc, jic+1 ) ! avoid out of bounds for tabin below 254 jjc1 = MIN( nyc, jjc+1 ) ! -- 255 256 tabout(ji,jj) = ( Bx * tabin(jic1,jjc ) + Ax * tabin(jic,jjc ) ) * Ay + & 257 & ( Bx * tabin(jic1,jjc1) + Ax * tabin(jic,jjc1) ) * By 252 258 END DO 253 259 END DO -
utils/tools_MERGE_2019/NESTING/src/agrif_partial_steps.f90
r10025 r12412 174 174 END DO 175 175 ! 176 !!$ IF( partial_steps ) THEN 176 177 ! Compute bathy_level for ocean points (i.e. the number of ocean levels) 177 178 ! find the number of ocean levels such that the last level thickness 178 179 ! is larger than the minimum of e3zps_min and e3zps_rat * e3t (where 179 180 ! e3t is the reference level thickness 180 !181 181 DO jk = N-1, 1, -1 182 182 zdepth = gdepw(jk) + MIN( e3zps_min, e3t(jk)*e3zps_rat ) … … 184 184 DO ji = 1, jpi 185 185 IF( 0. < Grid%bathy_meter(ji,jj) .AND. Grid%bathy_meter(ji,jj) <= zdepth ) & 186 186 Grid%bathy_level(ji,jj) = jk-1 187 187 END DO 188 188 END DO 189 189 END DO 190 190 !!$ ELSE 191 !!$ DO jj = 1,jpj 192 !!$ DO ji = 1,jpi 193 !!$ ! 194 !!$ IF (Grid%bathy_meter(ji,jj) .EQ. 0.0 ) THEN 195 !!$ Grid%bathy_level(ji,jj)=0 196 !!$ ELSE 197 !!$ ! 198 !!$ k1=2 ! clem: minimum levels = 2 ??? 199 !!$ DO WHILE (k1 .LT. (N-1)) 200 !!$ IF ((Grid%bathy_meter(ji,jj).GE.gdepw(k1)) & 201 !!$ .AND.(Grid%bathy_meter(ji,jj).LE.gdepw(k1+1))) EXIT 202 !!$ k1=k1+1 203 !!$ END DO 204 !!$ Grid%bathy_level(ji,jj)=k1 205 !!$ ! 206 !!$ ENDIF 207 !!$ ! 208 !!$ END DO 209 !!$ END DO 210 !!$ 211 !!$ ENDIF 191 212 192 213 CALL bathymetry_control(grid%bathy_level) … … 675 696 SUBROUTINE bathymetry_control(mbathy) 676 697 677 INTEGER :: i,j, jl698 INTEGER :: ji, jj, jl 678 699 INTEGER :: icompt, ibtest, ikmax 679 700 REAL*8, DIMENSION(:,:) :: mbathy … … 692 713 DO jl = 1, 2 693 714 ! 694 DO j = 2, SIZE(mbathy,2)-1695 DO i = 2, SIZE(mbathy,1)-1696 697 ibtest = MAX( mbathy( i-1,j), mbathy(i+1,j),mbathy(i,j-1),mbathy(i,j+1) )715 DO jj = 2, SIZE(mbathy,2)-1 716 DO ji = 2, SIZE(mbathy,1)-1 717 718 ibtest = MAX( mbathy(ji-1,jj), mbathy(ji+1,jj),mbathy(ji,jj-1),mbathy(ji,jj+1) ) 698 719 ! 699 IF( ibtest < mbathy( i,j) ) THEN720 IF( ibtest < mbathy(ji,jj) ) THEN 700 721 ! 701 WRITE(*,*) 'grid-point(i,j)= ', i,j,'is changed from',mbathy(i,j),' to ', ibtest702 mbathy( i,j) = ibtest722 WRITE(*,*) 'grid-point(i,j)= ',ji,jj,'is changed from',mbathy(ji,jj),' to ', ibtest 723 mbathy(ji,jj) = ibtest 703 724 icompt = icompt + 1 704 725 ! … … 720 741 721 742 ikmax = 0 722 DO j = 1, SIZE(mbathy,2)743 DO jj = 1, SIZE(mbathy,2) 723 744 DO ji = 1, SIZE(mbathy,1) 724 ikmax = MAX( ikmax, NINT(mbathy( i,j)) )745 ikmax = MAX( ikmax, NINT(mbathy(ji,jj)) ) 725 746 END DO 726 747 END DO -
utils/tools_MERGE_2019/NESTING/src/agrif_readwrite.f90
r11231 r12412 241 241 CALL write_ncdf_var('nav_lon' ,dimnames,name,Grid%nav_lon ,'float') 242 242 CALL write_ncdf_var('nav_lat' ,dimnames,name,Grid%nav_lat ,'float') 243 CALL write_ncdf_var(parent_level_name,dimnames,name, Grid%bathy_level,'float')243 CALL write_ncdf_var(parent_level_name,dimnames,name,NINT(Grid%bathy_level),'integer') 244 244 ! 245 245 CALL copy_ncdf_att('nav_lon' ,TRIM(parent_bathy_level),name,MINVAL(Grid%nav_lon),MAXVAL(Grid%nav_lon)) … … 545 545 ln_sco = 0 546 546 ln_isfcav = 0 547 IF( partial_steps ) THEN547 !!$ IF( partial_steps ) THEN 548 548 ln_zps = 1 549 549 ln_zco = 0 550 CALL zgr_zps( nx, ny, gdept_1d, gdepw_1d, e3t_1d, e3w_1d, zbathy, & ! input 551 & mbathy, e3t_0, e3u_0, e3v_0, e3f_0, e3w_0, e3uw_0, e3vw_0 ) ! output 552 ELSE 553 ln_zps = 0 554 ln_zco = 1 555 ENDIF 550 !!$ ELSE 551 !!$ ln_zps = 0 552 !!$ ln_zco = 1 553 !!$ ENDIF 554 555 CALL zgr_zps( nx, ny, gdept_1d, gdepw_1d, e3t_1d, e3w_1d, zbathy, & ! input 556 & mbathy, e3t_0, e3u_0, e3v_0, e3f_0, e3w_0, e3uw_0, e3vw_0 ) ! output 556 557 557 558 ! top/bottom levels … … 865 866 zmax = gdepw_1d(N) + e3t_1d(N) ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_1d(N-1) ) 866 867 zbathy(:,:) = MIN( zmax , zbathy(:,:) ) ! bounded value of bathy (min already set at the end of zgr_bat) 867 WHERE( zbathy(:,:) == 0. ) ; mbathy(:,:) = 0 868 WHERE( zbathy(:,:) == 0. ) ; mbathy(:,:) = 0 ! land : set mbathy to 0 868 869 ELSEWHERE ; mbathy(:,:) = N-1 ! ocean : initialize mbathy to the max ocean level 869 870 END WHERE 870 871 871 ! Compute mbathy for ocean points (i.e. the number of ocean levels) 872 ! find the number of ocean levels such that the last level thickness 873 ! is larger than the minimum of e3zps_min and e3zps_rat * e3t_1d (where 874 ! e3t_1d is the reference level thickness 875 DO jk = N-1, 1, -1 876 zdepth = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 877 WHERE( 0. < zbathy(:,:) .AND. zbathy(:,:) <= zdepth ) mbathy(:,:) = jk-1 878 END DO 879 872 !!$ IF( partial_steps ) THEN 873 ! Compute mbathy for ocean points (i.e. the number of ocean levels) 874 ! find the number of ocean levels such that the last level thickness 875 ! is larger than the minimum of e3zps_min and e3zps_rat * e3t_1d (where 876 ! e3t_1d is the reference level thickness 877 DO jk = N-1, 1, -1 878 zdepth = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 879 WHERE( 0. < zbathy(:,:) .AND. zbathy(:,:) <= zdepth ) mbathy(:,:) = jk-1 880 END DO 881 !!$ ELSE 882 !!$ DO jk = 1, N 883 !!$ WHERE( 0. < zbathy(:,:) .AND. zbathy(:,:) >= gdepw_1d(jk) ) mbathy(:,:) = jk-1 884 !!$ END DO 885 !!$ ENDIF 886 880 887 ! Scale factors and depth at T- and W-points 881 888 DO jk = 1, N ! intitialization to the reference z-coordinate … … 891 898 ! 892 899 ! Scale factors and depth at T- and W-points 893 ! IF ( ln_isfcav == 0 ) THEN 894 DO jj = 1, ny 895 DO ji = 1, nx 896 ik = mbathy(ji,jj) 897 IF( ik > 0 ) THEN ! ocean point only 898 ! max ocean level case 899 IF( ik == N-1 ) THEN 900 zdepwp = zbathy(ji,jj) 901 ze3tp = zbathy(ji,jj) - gdepw_1d(ik) 902 ze3wp = 0.5 * e3w_1d(ik) * ( 1. + ( ze3tp/e3t_1d(ik) ) ) 903 e3t_0(ji,jj,ik ) = ze3tp 904 e3t_0(ji,jj,ik+1) = ze3tp 905 e3w_0(ji,jj,ik ) = ze3wp 906 e3w_0(ji,jj,ik+1) = ze3tp 907 gdepw_0(ji,jj,ik+1) = zdepwp 908 gdept_0(ji,jj,ik ) = gdept_1d(ik-1) + ze3wp 909 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp 910 ! 911 ELSE ! standard case 912 IF( zbathy(ji,jj) <= gdepw_1d(ik+1) ) THEN ; gdepw_0(ji,jj,ik+1) = zbathy(ji,jj) 913 ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 914 ENDIF 915 !gm Bug? check the gdepw_1d 916 ! ... on ik 917 gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) & 918 & * ((gdept_1d( ik ) - gdepw_1d(ik) ) & 919 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) )) 920 e3t_0 (ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) & 921 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) ) 922 e3w_0(ji,jj,ik) = 0.5 * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2. * gdepw_1d(ik) ) & 923 & * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 924 ! ... on ik+1 925 e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 926 e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 927 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik) 900 ! IF ( ln_isfcav == 0 ) THEN 901 DO jj = 1, ny 902 DO ji = 1, nx 903 ik = mbathy(ji,jj) 904 IF( ik > 0 ) THEN ! ocean point only 905 ! max ocean level case 906 IF( ik == N-1 ) THEN 907 zdepwp = zbathy(ji,jj) 908 ze3tp = zbathy(ji,jj) - gdepw_1d(ik) 909 ze3wp = 0.5 * e3w_1d(ik) * ( 1. + ( ze3tp/e3t_1d(ik) ) ) 910 e3t_0(ji,jj,ik ) = ze3tp 911 e3t_0(ji,jj,ik+1) = ze3tp 912 e3w_0(ji,jj,ik ) = ze3wp 913 e3w_0(ji,jj,ik+1) = ze3tp 914 gdepw_0(ji,jj,ik+1) = zdepwp 915 gdept_0(ji,jj,ik ) = gdept_1d(ik-1) + ze3wp 916 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp 917 ! 918 ELSE ! standard case 919 IF( zbathy(ji,jj) <= gdepw_1d(ik+1) ) THEN ; gdepw_0(ji,jj,ik+1) = zbathy(ji,jj) 920 ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 928 921 ENDIF 922 !gm Bug? check the gdepw_1d 923 ! ... on ik 924 gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) & 925 & * ((gdept_1d( ik ) - gdepw_1d(ik) ) & 926 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) )) 927 e3t_0 (ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) & 928 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) ) 929 e3w_0(ji,jj,ik) = 0.5 * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2. * gdepw_1d(ik) ) & 930 & * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 931 ! ... on ik+1 932 e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 933 e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 934 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik) 929 935 ENDIF 930 END DO936 ENDIF 931 937 END DO 932 ! 933 ! DO jj = 1, ny 934 ! DO ji = 1, nx 935 ! ik = mbathy(ji,jj) 936 ! IF( ik > 0 ) THEN ! ocean point only 937 ! e3tp (ji,jj) = e3t_0(ji,jj,ik) 938 ! e3wp (ji,jj) = e3w_0(ji,jj,ik) 939 ! ENDIF 940 ! END DO 941 ! END DO 942 ! END IF 938 END DO 939 ! 940 ! DO jj = 1, ny 941 ! DO ji = 1, nx 942 ! ik = mbathy(ji,jj) 943 ! IF( ik > 0 ) THEN ! ocean point only 944 ! e3tp (ji,jj) = e3t_0(ji,jj,ik) 945 ! e3wp (ji,jj) = e3w_0(ji,jj,ik) 946 ! ENDIF 947 ! END DO 948 ! END DO 949 ! END IF 943 950 ! 944 951 ! Scale factors and depth at U-, V-, UW and VW-points -
utils/tools_MERGE_2019/NESTING/src/agrif_types.f90
r11206 r12412 192 192 ENDIF 193 193 ! 194 IF( .NOT.partial_steps ) THEN 195 WRITE(*,*) 'E R R O R: partial_steps must be set to true otherwise very thin bottom layers can be created' 196 STOP 197 ENDIF 198 ! 194 199 ELSE 195 200 ! -
utils/tools_MERGE_2019/SIREN/src/create_boundary.F90
r12080 r12412 937 937 ELSE 938 938 939 WRITE(cl_errormsg,*) " ERROR : can't find "//TRIM(c l_namelist)939 WRITE(cl_errormsg,*) " ERROR : can't find "//TRIM(cd_namelist) 940 940 CALL fct_help(cp_myname,cl_errormsg) 941 941 CALL EXIT(1) -
utils/tools_MERGE_2019/SIREN/src/docsrc/5_changeLog.md
r12080 r12412 5 5 # Release 2019-12-03 {#rev_2019-12-03} 6 6 ## New features 7 - M Siren/src/iom_cdf.f90 :7 - M src/iom_cdf.f90 : 8 8 - write netcdf file as netcdf4 9 9 10 10 # Release 2019-11-05 {#rev_2019-11-05} 11 11 ## New features 12 - M Siren/src/function.f9012 - M src/function.f90 13 13 - M src/create_bathy.f90 : 14 14 - add help and version optional arguments
Note: See TracChangeset
for help on using the changeset viewer.
