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 12412 for utils/tools_MERGE_2019/NESTING/src/agrif_create_bathy.f90 – NEMO

Ignore:
Timestamp:
2020-02-19T18:48:10+01:00 (4 years ago)
Author:
smueller
Message:

Synchronisation of the 2019 development branch for the tools directory (/utils/tools) with /utils/tools@12273

File:
1 edited

Legend:

Unmodified
Added
Removed

  • utils/tools_MERGE_2019/NESTING/src/agrif_create_bathy.f90

    r10381 r12412  
    8484     STOP 
    8585  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 
    97111        status = Read_bathy_level(TRIM(parent_bathy_level),G0) 
    98112        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      
    128114     ! change longitudes (from -180:180 to 0:360) 
    129      IF ( Pacifique ) THEN 
     115     IF(Pacifique) THEN 
    130116        WHERE(G0%nav_lon < 0.001)   G0%nav_lon = G0%nav_lon + 360. 
    131117     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 
    279164     !                                                               ! -----------------------------  
    280165     IF( type_bathy_interp == 0 .OR. type_bathy_interp == 1 ) THEN   ! arithmetic or median averages 
    … …  
    375260           type_bathy_interp = 2 
    376261        ENDIF 
    377             
     262 
     263        DEALLOCATE(trouble_points) 
     264 
    378265     ENDIF 
    379266     !                                                       ! -----------------------------  
    380267     IF( type_bathy_interp == 2) THEN                        ! Bilinear interpolation 
    381268        !                                                    ! -----------------------------  
    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 
    415276        !             
    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 
    420359     ! 
    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 
    423375     ! 
    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. 
    434379     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 
    493409     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 
    509417     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)) 
    536424        END DO 
    537  
    538         wghts = 1. 
    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 
     425     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)) 
    544432        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. 
    603475                 ENDIF 
    604               ENDDO 
     476 
     477              ENDIF 
    605478           ENDDO 
    606479        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 
    631501  ! 
    632502END PROGRAM create_bathy 
Note: See TracChangeset for help on using the changeset viewer.