New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 12264 for utils/tools/NESTING/src/agrif_create_bathy.f90 – NEMO

Ignore:
Timestamp:
2019-12-17T18:32:39+01:00 (4 years ago)
Author:
clem
Message:

finishing cleaning the nesting tools (as far as I could). I will stop there

File:
1 edited

Legend:

Unmodified
Added
Removed
  • utils/tools/NESTING/src/agrif_create_bathy.f90

    r12259 r12264  
    105105     DEALLOCATE( G0%nav_lon, G0%nav_lat ) 
    106106     status = read_bathy_coord(TRIM(elevation_database),G0,G1,Pacifique) 
    107      CALL meter_to_levels(G0) 
    108107  ELSE                                                      ! read G0%bathy_meter (on the global grid) 
    109      IF( TRIM(parent_bathy_level) /= '' ) THEN 
     108     IF( TRIM(parent_bathy_meter) /= '') THEN 
     109        status = read_bathy_meter(TRIM(parent_bathy_meter),G0) 
     110     ELSE 
    110111        status = Read_bathy_level(TRIM(parent_bathy_level),G0) 
    111112        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 
     113     ENDIF      
    116114     ! change longitudes (from -180:180 to 0:360) 
    117115     IF(Pacifique) THEN 
     
    139137  ENDIF 
    140138   
    141   !                                                   ! ------------------------------------------------------------------ 
    142   IF( .NOT.new_topo .AND. .NOT.partial_steps ) THEN   ! **** First option : no new topo file & no partial steps 
    143      !                                                !                 ==> interpolate levels directly from parent file 
    144      !                                                !                     using bilinear interpolation 
    145      !                                                ! ------------------------------------------------------------------ 
    146      WRITE(*,*) '*** First option: no new topo file & no partial steps' 
    147      ! 
    148      IF( .NOT. identical_grids ) THEN 
    149         !            
    150         WRITE(*,*) 'Interpolation of high resolution bathymetry on child grid: Bilinear interpolation ...' 
     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 
     164     !                                                               ! -----------------------------  
     165     IF( type_bathy_interp == 0 .OR. type_bathy_interp == 1 ) THEN   ! arithmetic or median averages 
     166        !                                                            ! -----------------------------  
     167        ALLOCATE(trouble_points(nxfin,nyfin)) 
     168        trouble_points(:,:) = 0 
     169        !                        
     170        DO jj = 2, nyfin 
     171           DO ji = 2, nxfin 
     172              !        
     173              ! fine grid cell extension                
     174              Cell_lonmin = MIN(G1%glamf(ji-1,jj-1),G1%glamf(ji,jj-1),G1%glamf(ji,jj),G1%glamf(ji-1,jj)) 
     175              Cell_lonmax = MAX(G1%glamf(ji-1,jj-1),G1%glamf(ji,jj-1),G1%glamf(ji,jj),G1%glamf(ji-1,jj)) 
     176              Cell_latmin = MIN(G1%gphif(ji-1,jj-1),G1%gphif(ji,jj-1),G1%gphif(ji,jj),G1%gphif(ji-1,jj)) 
     177              Cell_latmax = MAX(G1%gphif(ji-1,jj-1),G1%gphif(ji,jj-1),G1%gphif(ji,jj),G1%gphif(ji-1,jj))  
     178              !                
     179              ! look for points in G0 (bathy dataset) contained in the fine grid cells   
     180              iimin = 1 
     181              DO WHILE( G0%nav_lon(iimin,1) < Cell_lonmin )  
     182                 iimin = iimin + 1 
     183              ENDDO 
     184              !                
     185              jjmin = 1 
     186              DO WHILE( G0%nav_lat(iimin,jjmin) < Cell_latmin )  
     187                 jjmin = jjmin + 1 
     188              ENDDO 
     189              !                 
     190              iimax = iimin  
     191              DO WHILE( G0%nav_lon(iimax,1) <= Cell_lonmax )  
     192                 iimax = iimax + 1 
     193                 iimax = MIN( iimax,SIZE(G0%bathy_meter,1)) 
     194              ENDDO 
     195              !                                
     196              jjmax = jjmin  
     197              DO WHILE( G0%nav_lat(iimax,jjmax) <= Cell_latmax )  
     198                 jjmax = jjmax + 1 
     199                 jjmax = MIN( jjmax,SIZE(G0%bathy_meter,2)) 
     200              ENDDO 
     201              ! 
     202              IF( ln_agrif_domain ) THEN 
     203                 iimax = iimax-1 
     204                 jjmax = jjmax-1 
     205              ELSE 
     206                 iimax = MAX(iimin,iimax-1) 
     207                 jjmax = MAX(jjmin,jjmax-1) 
     208              ENDIF 
     209              !                
     210              iimin = MAX( iimin,1 ) 
     211              jjmin = MAX( jjmin,1 ) 
     212              iimax = MIN( iimax,SIZE(G0%bathy_meter,1)) 
     213              jjmax = MIN( jjmax,SIZE(G0%bathy_meter,2)) 
     214 
     215              nxhr = iimax - iimin + 1 
     216              nyhr = jjmax - jjmin + 1                     
     217 
     218              IF( nxhr == 0 .OR. nyhr == 0 ) THEN 
     219                 ! 
     220                 trouble_points(ji,jj) = 1 
     221                 ! 
     222              ELSE 
     223                 ! 
     224                 ALLOCATE( vardep(nxhr,nyhr), mask_oce(nxhr,nyhr) ) 
     225                 vardep(:,:) = G0%bathy_meter(iimin:iimax,jjmin:jjmax) 
     226                 ! 
     227                 WHERE( vardep(:,:) .GT. 0. )   ;   mask_oce = 1 ; 
     228                 ELSEWHERE                      ;   mask_oce = 0 ; 
     229                 ENDWHERE 
     230                 ! 
     231                 nxyhr = nxhr*nyhr 
     232                 IF( SUM(mask_oce) < 0.5*(nxyhr) ) THEN ! if more than half of the points are on land then bathy fine = 0 
     233                    G1%bathy_meter(ji,jj) = 0. 
     234                 ELSE 
     235                    IF( type_bathy_interp == 0 ) THEN     ! Arithmetic average 
     236                       G1%bathy_meter(ji,jj) = SUM( vardep(:,:) * mask_oce(:,:) ) / SUM( mask_oce(:,:) ) 
     237                    ELSE                                  ! Median average 
     238                       ALLOCATE(vardep1d(nxyhr)) 
     239                       vardep1d = RESHAPE(vardep,(/ nxyhr /) ) 
     240                       !!CALL ssort(vardep1d,nxyhr) 
     241                       CALL quicksort(vardep1d,1,nxyhr) 
     242                       ! 
     243                       ! Calculate median 
     244                       IF (MOD(nxyhr,2) .NE. 0) THEN 
     245                          G1%bathy_meter(ji,jj) = vardep1d( nxyhr/2 + 1 ) 
     246                       ELSE 
     247                          G1%bathy_meter(ji,jj) = 0.5 * ( vardep1d(nxyhr/2) + vardep1d(nxyhr/2+1) ) 
     248                       END IF 
     249                       DEALLOCATE(vardep1d)    
     250                    ENDIF 
     251                 ENDIF 
     252                 DEALLOCATE (mask_oce,vardep) 
     253                 ! 
     254              ENDIF 
     255           ENDDO 
     256        ENDDO 
     257 
     258        IF( SUM( trouble_points ) > 0 ) THEN 
     259           PRINT*,'too much empty cells, proceed to bilinear interpolation' 
     260           type_bathy_interp = 2 
     261        ENDIF 
     262 
     263        DEALLOCATE(trouble_points) 
     264 
     265     ENDIF 
     266     !                                                       ! -----------------------------  
     267     IF( type_bathy_interp == 2) THEN                        ! Bilinear interpolation 
     268        !                                                    ! -----------------------------  
     269 
     270        ALLOCATE(masksrc(SIZE(G0%bathy_meter,1),SIZE(G0%bathy_meter,2))) 
     271        ALLOCATE(bathy_test(nxfin,nyfin)) 
    151272        ! 
    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.   ; 
     273        WHERE(G0%bathy_meter.LE.0)   ;   masksrc = .FALSE.   ; 
    155274        ELSEWHERE                    ;   masksrc = .TRUE.    ; 
    156275        END WHERE 
    157         ! 
    158         ! compute remapping matrix thanks to SCRIP package 
     276        !             
     277        ! compute remapping matrix thanks to SCRIP package             
    159278        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)   
     279        CALL make_remap(G0%bathy_meter,bathy_test,nxfin,nyfin,matrix,src_add,dst_add)   
     280        !                                   
     281        G1%bathy_meter = bathy_test                
     282        !             
    161283        DEALLOCATE(masksrc) 
    162         ! 
    163      ENDIF 
    164       
    165      ! compute constant bathymetry for Parent-Child bathymetry connection 
    166      CALL init_constant_bathy(G0%bathy_meter,bathy_fin_constant) 
    167      ! 
    168      ! replace child bathymetry by parent bathymetry at the boundaries 
    169      IF( ln_agrif_domain ) THEN 
    170         boundary = npt_copy*irafx + nbghostcellsfine + 1  
    171      ELSE 
    172         boundary = npt_copy*irafx 
    173      ENDIF 
    174      G1%bathy_meter(1:boundary,:) = bathy_fin_constant(1:boundary,:) 
    175      G1%bathy_meter(:,1:boundary) = bathy_fin_constant(:,1:boundary) 
    176      G1%bathy_meter(nxfin-boundary+1:nxfin,:) = bathy_fin_constant(nxfin-boundary+1:nxfin,:) 
    177      G1%bathy_meter(:,nyfin-boundary+1:nyfin) = bathy_fin_constant(:,nyfin-boundary+1:nyfin) 
    178      DEALLOCATE(bathy_fin_constant) 
    179      !                   
    180      ! bathymetry smoothing (everywhere except at the boundaries)  
    181      IF( smoothing ) THEN 
    182         IF( ln_agrif_domain ) THEN 
    183            CALL smooth_topo(G1%bathy_meter(boundary:nxfin-boundary+1,boundary:nyfin-boundary+1),nbiter) 
    184         ELSE 
    185            CALL smooth_topo(G1%bathy_meter(boundary+1:nxfin-boundary,boundary+1:nyfin-boundary),nbiter) 
    186         ENDIF 
    187      ENDIF 
    188      ! 
    189      ! From meters to levels 
    190      CALL meter_to_levels(G1) 
    191      G1%bathy_level=NINT(G1%bathy_level) 
    192      !        
    193      ! Check errors in bathy 
    194      DO jj=1,nyfin 
    195         DO ji=1,nxfin 
    196            IF (G1%bathy_level(ji,jj).LT.0) THEN 
    197               PRINT *,'error in ',ji,jj,G1%bathy_level(ji,jj) 
    198               STOP 
    199            ENDIF 
    200         ENDDO 
    201      ENDDO 
    202      WHERE ((G1%bathy_level.LT.3).AND.(G1%bathy_level.GT.0))   G1%bathy_level=3 
    203      ! 
    204      ! remove closed seas 
    205      IF (removeclosedseas) THEN 
    206         ALLOCATE(bathy_test(nxfin,nyfin)) 
    207  
    208         bathy_test=0 
    209         WHERE (G1%bathy_level(1,:)     .GT.0)   bathy_test(1,:)=1 
    210         WHERE (G1%bathy_level(nxfin,:) .GT.0)   bathy_test(nxfin,:)=1 
    211         WHERE (G1%bathy_level(:,1)     .GT.0)   bathy_test(:,1)=1 
    212         WHERE (G1%bathy_level(:,nyfin) .GT.0)   bathy_test(:,nyfin)=1 
    213  
    214         nbadd = 1 
    215         DO WHILE (nbadd.NE.0) 
    216            nbadd = 0 
    217            DO jj=2,nyfin-1 
    218               DO ji=2,nxfin-1 
    219                  IF (G1%bathy_level(ji,jj).GT.0) THEN 
    220                     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 
    221                        IF (bathy_test(ji,jj).NE.1) nbadd = nbadd + 1 
    222                        bathy_test(ji,jj)=1 
    223                     ENDIF 
    224                  ENDIF 
    225               ENDDO 
    226            ENDDO 
    227         ENDDO 
    228  
    229         WHERE (bathy_test.EQ.0.)   G1%bathy_level = 0 
    230  
    231         DEALLOCATE(bathy_test) 
    232      ENDIF 
    233      ! 
    234      CALL levels_to_meter(G1) ! needed for domcfg 
    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      ! 
    243      ! store interpolation result in output file 
    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) 
    247      !!status = write_domcfg(TRIM(parent_domcfg_out),G0) ! do not activate it 
    248  
    249      WRITE(*,*) '****** Bathymetry successfully created if no new topo and no partial steps ******' 
    250      STOP 
    251      ! 
    252      !                                                ! ----------------------------------------------------- 
    253   ELSE                                                ! **** Second option : new topo file or partial steps      
    254      !                                                ! -----------------------------------------------------  
    255      WRITE(*,*) '*** Second option : new topo or partial steps' 
    256      ! 
    257      ! === Here: G0 is the grid associated with the new topography (as gebco or etopo) === 
    258      ! 
    259      ! what type of interpolation for bathymetry 
    260      IF( type_bathy_interp == 0 ) THEN 
    261         WRITE(*,*) 'Interpolation of high resolution bathymetry on child grid: Arithmetic average ...' 
    262      ELSE IF( type_bathy_interp == 1 ) THEN 
    263         WRITE(*,*) 'Interpolation of high resolution bathymetry on child grid: Median average ...' 
    264      ELSE IF( type_bathy_interp == 2 ) THEN      
    265         WRITE(*,*) 'Interpolation of high resolution bathymetry on child grid: Bilinear interpolation ...' 
    266      ELSE      
    267         WRITE(*,*) 'bad value for type_bathy_interp variable ( must be 0, 1 or 2 )' 
    268         STOP  
    269      ENDIF 
    270      ! 
    271      ! 
    272      ! --------------------------------------------------------------------------------- 
    273      ! ===                 Bathymetry of the fine grid (step1)                       === 
    274      ! --------------------------------------------------------------------------------- 
    275      ! ==> It gives G1%bathy_meter from G0%bathy_meter 
    276      ! --------------------------------------------------------------------------------- 
    277  
    278      IF( .NOT. identical_grids ) THEN 
    279         !                                                               ! -----------------------------  
    280         IF( type_bathy_interp == 0 .OR. type_bathy_interp == 1 ) THEN   ! arithmetic or median averages 
    281            !                                                            ! -----------------------------  
    282            ALLOCATE(trouble_points(nxfin,nyfin)) 
    283            trouble_points(:,:) = 0 
    284            !                        
    285            DO jj = 2, nyfin 
    286               DO ji = 2, nxfin 
    287                  !        
    288                  ! fine grid cell extension                
    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)) 
    292                  Cell_latmax = MAX(G1%gphif(ji-1,jj-1),G1%gphif(ji,jj-1),G1%gphif(ji,jj),G1%gphif(ji-1,jj))  
    293                  !                
    294                  ! look for points in G0 (bathy dataset) contained in the fine grid cells   
    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 
    308                     iimax = MIN( iimax,SIZE(G0%bathy_meter,1)) 
    309                  ENDDO 
    310                  !                                
    311                  jjmax = jjmin  
    312                  DO WHILE( G0%nav_lat(iimax,jjmax) <= Cell_latmax )  
    313                     jjmax = jjmax + 1 
    314                     jjmax = MIN( jjmax,SIZE(G0%bathy_meter,2)) 
    315                  ENDDO 
    316                  ! 
    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 
    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 
    334                     ! 
    335                     trouble_points(ji,jj) = 1 
    336                     ! 
    337                  ELSE 
    338                     ! 
    339                     ALLOCATE( vardep(nxhr,nyhr), mask_oce(nxhr,nyhr) ) 
    340                     vardep(:,:) = G0%bathy_meter(iimin:iimax,jjmin:jjmax) 
    341                     ! 
    342                     WHERE( vardep(:,:) .GT. 0. )   ;   mask_oce = 1 ; 
    343                     ELSEWHERE                      ;   mask_oce = 0 ; 
    344                     ENDWHERE 
    345                     ! 
    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 
    348                        G1%bathy_meter(ji,jj) = 0. 
    349                     ELSE 
    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 
    353                           ALLOCATE(vardep1d(nxyhr)) 
    354                           vardep1d = RESHAPE(vardep,(/ nxyhr /) ) 
    355                           !!CALL ssort(vardep1d,nxyhr) 
    356                           CALL quicksort(vardep1d,1,nxyhr) 
    357                           ! 
    358                           ! Calculate median 
    359                           IF (MOD(nxyhr,2) .NE. 0) THEN 
    360                              G1%bathy_meter(ji,jj) = vardep1d( nxyhr/2 + 1 ) 
    361                           ELSE 
    362                              G1%bathy_meter(ji,jj) = 0.5 * ( vardep1d(nxyhr/2) + vardep1d(nxyhr/2+1) ) 
    363                           END IF 
    364                           DEALLOCATE(vardep1d)    
    365                        ENDIF 
    366                     ENDIF 
    367                     DEALLOCATE (mask_oce,vardep) 
    368                     ! 
    369                  ENDIF 
    370               ENDDO 
    371            ENDDO 
    372  
    373            IF( SUM( trouble_points ) > 0 ) THEN 
    374               PRINT*,'too much empty cells, proceed to bilinear interpolation' 
    375               type_bathy_interp = 2 
    376            ENDIF 
    377  
    378            DEALLOCATE(trouble_points) 
    379  
    380         ENDIF 
    381         !                                                       ! -----------------------------  
    382         IF( type_bathy_interp == 2) THEN                        ! Bilinear interpolation 
    383            !                                                    ! -----------------------------  
    384  
    385            ALLOCATE(masksrc(SIZE(G0%bathy_meter,1),SIZE(G0%bathy_meter,2))) 
    386            ALLOCATE(bathy_test(nxfin,nyfin)) 
    387            ! 
    388            WHERE(G0%bathy_level.LE.0)   ;   masksrc = .FALSE.   ; 
    389            ELSEWHERE                    ;   masksrc = .TRUE.    ; 
    390            END WHERE 
    391            !             
    392            ! compute remapping matrix thanks to SCRIP package             
    393            CALL get_remap_matrix(G0%nav_lat,G1%nav_lat,G0%nav_lon,G1%nav_lon,masksrc,matrix,src_add,dst_add) 
    394            CALL make_remap(G0%bathy_meter,bathy_test,nxfin,nyfin,matrix,src_add,dst_add)   
    395            !                                   
    396            G1%bathy_meter = bathy_test                
    397            !             
    398            DEALLOCATE(masksrc) 
    399            DEALLOCATE(bathy_test)  
    400             
    401         ENDIF 
    402         !             
    403      ENDIF ! not identical grids 
    404      ! --- 
    405      ! At this stage bathymetry in meters has already been interpolated on fine grid 
    406      !                    => G1%bathy_meter(nxfin,nyfin) 
    407      ! 
    408      ! Also G0 was the grid from the new bathymetry data set (etopo, gebco...) and not the coarse grid 
    409      ! --- 
    410      ! 
    411      ! --------------------------------------------------------------------------------- 
    412      ! ===                 Bathymetry of the fine grid (step2)                       === 
    413      ! --------------------------------------------------------------------------------- 
    414      ! ==> It gives an update of G1%bathy_meter and G1%bathy_level 
    415      ! --------------------------------------------------------------------------------- 
    416      ! From here on: G0 is the coarse grid 
    417      ! 
    418      ! Coarse grid bathymetry : G0%bathy_meter (on the global grid) 
    419      IF( TRIM(parent_bathy_meter) /= '') THEN 
    420         status = read_bathy_meter(TRIM(parent_bathy_meter),G0) 
    421      ELSE 
    422         status = Read_bathy_level(TRIM(parent_bathy_level),G0) 
    423         CALL levels_to_meter(G0) 
    424      ENDIF 
    425       
    426      ! Coarse grid coordinatees : G0 coordinates 
    427      DEALLOCATE(G0%nav_lat,G0%nav_lon) 
    428      status = Read_coordinates(TRIM(parent_coordinate_file),G0) 
    429  
    430      ! allocate temporary arrays                   
    431      IF (.NOT.ASSOCIATED(G0%gdepw_ps))       ALLOCATE(G0%gdepw_ps    (SIZE(G0%bathy_meter,1),SIZE(G0%bathy_meter,2))) 
    432      IF (.NOT.ASSOCIATED(G1%gdepw_ps))       ALLOCATE(G1%gdepw_ps    (SIZE(G1%bathy_meter,1),SIZE(G1%bathy_meter,2))) 
    433      IF (.NOT.ASSOCIATED(gdepw_ps_interp))   ALLOCATE(gdepw_ps_interp(SIZE(G1%bathy_meter,1),SIZE(G1%bathy_meter,2))) 
    434      !                        
    435      IF( ln_agrif_domain ) THEN 
    436         boundary = npt_copy*irafx + nbghostcellsfine + 1 
    437      ELSE 
    438         boundary = npt_copy*irafx 
    439      ENDIF 
    440      ! 
    441      ! compute G0%gdepw_ps and G1%gdepw_ps 
    442      CALL get_partial_steps(G0)  
    443      CALL get_partial_steps(G1) 
    444      CALL bathymetry_control(G0%Bathy_level) 
    445  
    446      ! --------------------------------------- 
    447      ! Bathymetry at the boundaries (npt_copy)                       
    448      ! --------------------------------------- 
    449      ! 1st step: interpolate coarse bathy on the fine grid (using partial steps or not) 
    450      IF( partial_steps .AND. ln_agrif_domain ) THEN                    
    451         CALL Check_interp(G0,gdepw_ps_interp) 
    452      ELSE 
    453         gdepw_ps_interp = 0. * G1%gdepw_ps 
    454         !!CALL agrif_interp(G0%gdepw_ps,gdepw_ps_interp,'T') 
    455         CALL init_constant_bathy(G0%gdepw_ps,gdepw_ps_interp) 
    456      ENDIF 
    457  
    458      IF (.NOT.ASSOCIATED(G1%wgt))   ALLOCATE(G1%wgt(SIZE(G1%bathy_meter,1),SIZE(G1%bathy_meter,2))) 
    459      G1%wgt(:,:) = 0. 
    460      IF ((.NOT.ASSOCIATED(G0%wgt)).AND.bathy_update) THEN  
    461         ALLOCATE(G0%wgt(SIZE(G0%nav_lat,1),SIZE(G0%nav_lat,2))) 
    462         G0%wgt(:,:) = 0. 
    463      ENDIF 
    464      ! 
     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 
    465351     ! 2nd step: copy parent bathymetry at the boundaries 
    466352     DO jj=1,nyfin   ! West and East 
     
    496382     END DO 
    497383     ! 
    498      !clem: recalculate interpolation everywhere before linear connection (useless to me) 
    499      IF( partial_steps .AND. ln_agrif_domain ) THEN                 
     384     !clem: recalculate interpolation everywhere before linear connection (useless to me??) 
     385     IF( ln_agrif_domain ) THEN                 
    500386        gdepw_ps_interp = 0. 
    501387        CALL Check_interp(G0,gdepw_ps_interp) 
     
    551437 
    552438     ENDIF 
    553  
    554      ! replace G1%bathy_meter by G1%gdepw_ps 
    555      G1%bathy_meter = G1%gdepw_ps 
    556      !                      
    557      ! -------------------- 
    558      ! Bathymetry smoothing                  
    559      ! -------------------- 
    560      IF( smoothing .AND. (.NOT.identical_grids) ) THEN 
    561         ! Chanut: smoothing everywhere then discard result in connection zone 
    562         CALL smooth_topo(G1%gdepw_ps(:,:),nbiter) 
    563         WHERE (G1%wgt(:,:)==0.) G1%bathy_meter(:,:) = G1%gdepw_ps(:,:) 
    564      ELSE 
    565         WRITE(*,*) 'No smoothing process only connection is carried out' 
    566      ENDIF 
    567      ! 
    568      ! ------------------ 
    569      ! Remove closed seas 
    570      ! ------------------ 
    571      IF (removeclosedseas) THEN 
    572         ALLOCATE(bathy_test(nxfin,nyfin)) 
    573         bathy_test=0. 
    574         WHERE (G1%bathy_meter(1,:)    .GT.0.)   bathy_test(1,:)=1 
    575         WHERE (G1%bathy_meter(nxfin,:).GT.0.)   bathy_test(nxfin,:)=1 
    576         WHERE (G1%bathy_meter(:,1)    .GT.0.)   bathy_test(:,1)=1 
    577         WHERE (G1%bathy_meter(:,nyfin).GT.0.)   bathy_test(:,nyfin)=1 
    578  
    579         nbadd = 1 
    580         DO WHILE (nbadd.NE.0) 
    581            nbadd = 0 
    582            DO jj=2,nyfin-1 
    583               DO ji=2,nxfin-1 
    584                  IF (G1%bathy_meter(ji,jj).GT.0.) THEN 
    585                     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 
    586                        IF (bathy_test(ji,jj).NE.1.) nbadd = nbadd + 1 
    587                        bathy_test(ji,jj)=1. 
    588                     ENDIF 
    589  
     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. 
    590475                 ENDIF 
    591               ENDDO 
     476 
     477              ENDIF 
    592478           ENDDO 
    593479        ENDDO 
    594         WHERE (bathy_test.EQ.0.)   G1%bathy_meter = 0. 
    595         DEALLOCATE(bathy_test) 
    596      ENDIF 
    597      ! 
    598      IF( partial_steps ) THEN 
    599         CALL get_partial_steps(G1)  ! recompute bathy_level and gdepw_ps for G1 (and correct bathy_meter) 
    600      ELSE 
    601         CALL meter_to_levels(G1)    ! convert bathymetry from meters to levels 
    602      ENDIF 
    603      ! 
    604      ! update parent grid 
    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      ! 
    611      ! store interpolation result in output file 
    612      IF( TRIM(parent_bathy_level) /= '' )   status = Write_Bathy_level(TRIM(child_level),G1) 
    613      IF( TRIM(parent_bathy_meter) /= '' )   status = Write_Bathy_meter(TRIM(child_meter),G1) 
    614      IF( TRIM(parent_domcfg_out)  /= '' )   status = write_domcfg(TRIM(child_domcfg),G1) 
    615      ! 
    616      WRITE(*,*) '****** Bathymetry successfully created ******' 
    617      STOP 
    618   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 
    619501  ! 
    620502END PROGRAM create_bathy 
Note: See TracChangeset for help on using the changeset viewer.