Changeset 12259


Ignore:
Timestamp:
2019-12-16T18:42:48+01:00 (10 months ago)
Author:
clem
Message:

partially clean the mess in the creation of bathymetry for agrif zoom. Last step in a next commit

File:
1 edited

Legend:

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

    r12253 r12259  
    8484     STOP 
    8585  ENDIF 
     86  ! 
     87  ! read fine and coarse grids coordinates file 
     88  status = Read_Coordinates(TRIM(parent_coordinate_file),G0) 
     89  status = Read_Coordinates(TRIM(child_coordinates),G1,Pacifique) 
     90  ! 
     91  ! check error in size 
     92  IF( imax > SIZE(G0%nav_lon,1) .OR. jmax > SIZE(G0%nav_lon,2) .OR. imax <= imin .OR. jmax <= jmin ) THEN 
     93     WRITE(*,*) 'ERROR ***** bad child grid definition ...' 
     94     WRITE(*,*) 'please check imin,jmin,imax,jmax,jpizoom,jpjzoom values'        
     95     STOP 
     96  ENDIF 
     97  IF( SIZE(G1%nav_lon,1) .NE. nxfin .OR. SIZE(G1%nav_lon,2) .NE. nyfin ) THEN 
     98     WRITE(*,*) 'ERROR ***** bad child coordinates file ...' 
     99     WRITE(*,*) 'please check that child coordinates file has been created with the same namelist' 
     100     STOP 
     101  ENDIF 
     102  ! 
     103  ! read bathymetry data set => G0%bathy_meter 
     104  IF( new_topo ) THEN                                       ! read G0%bathy_meter (on a reduced grid) and G1 coordinates 
     105     DEALLOCATE( G0%nav_lon, G0%nav_lat ) 
     106     status = read_bathy_coord(TRIM(elevation_database),G0,G1,Pacifique) 
     107     CALL meter_to_levels(G0) 
     108  ELSE                                                      ! read G0%bathy_meter (on the global grid) 
     109     IF( TRIM(parent_bathy_level) /= '' ) THEN 
     110        status = Read_bathy_level(TRIM(parent_bathy_level),G0) 
     111        CALL levels_to_meter(G0) 
     112     ELSEIF( TRIM(parent_bathy_level) == '' .AND. TRIM(parent_bathy_meter) /= '') THEN 
     113        status = read_bathy_meter(TRIM(parent_bathy_meter),G0) 
     114        CALL meter_to_levels(G0) 
     115     ENDIF 
     116     ! change longitudes (from -180:180 to 0:360) 
     117     IF(Pacifique) THEN 
     118        WHERE(G0%nav_lon < 0.001)   G0%nav_lon = G0%nav_lon + 360. 
     119     ENDIF 
     120  ENDIF 
     121  ! 
     122  ! 1st allocation of child grid bathy 
     123  ALLOCATE(G1%bathy_meter(nxfin,nyfin)) 
     124  G1%bathy_meter(:,:)=0. 
     125 
     126  ! check grids: if identical then do not interpolate 
     127  identical_grids = .FALSE. 
     128   
     129  IF(  SIZE(G0%nav_lat,1) == SIZE(G1%nav_lat,1) .AND. SIZE(G0%nav_lat,2) == SIZE(G1%nav_lat,2) .AND. & 
     130     & SIZE(G0%nav_lon,1) == SIZE(G1%nav_lon,1) .AND. SIZE(G0%nav_lon,2) == SIZE(G1%nav_lon,2) ) THEN 
     131     IF(  MAXVAL( ABS(G0%nav_lat(:,:)- G1%nav_lat(:,:)) ) < 0.0001 .AND. & 
     132        & MAXVAL( ABS(G0%nav_lon(:,:)- G1%nav_lon(:,:)) ) < 0.0001 ) THEN 
     133        WRITE(*,*) '' 
     134        WRITE(*,*) 'same grid between parent and child domains => NO INTERPOLATION' 
     135        WRITE(*,*) '' 
     136        G1%bathy_meter = G0%bathy_meter 
     137        identical_grids = .TRUE. 
     138     ENDIF 
     139  ENDIF 
     140   
    86141  !                                                   ! ------------------------------------------------------------------ 
    87142  IF( .NOT.new_topo .AND. .NOT.partial_steps ) THEN   ! **** First option : no new topo file & no partial steps 
     
    90145     !                                                ! ------------------------------------------------------------------ 
    91146     WRITE(*,*) '*** First option: no new topo file & no partial steps' 
    92      !       
    93      ! read coarse grid coordinates 
    94      status = Read_Coordinates(TRIM(parent_coordinate_file),G0)     
    95  
    96      ! read coarse grid bathymetry 
    97      IF( TRIM(parent_bathy_level) /= '' ) THEN 
    98         status = Read_bathy_level(TRIM(parent_bathy_level),G0) 
    99         CALL levels_to_meter(G0) 
    100      ELSEIF( TRIM(parent_bathy_level) == '' .AND. TRIM(parent_bathy_meter) /= '') THEN 
    101         status = read_bathy_meter(TRIM(parent_bathy_meter),G0) 
    102         CALL meter_to_levels(G0) 
    103      ENDIF 
    104      !            
    105      ! read fine grid coordinates 
    106      status = Read_Coordinates(TRIM(child_coordinates),G1,pacifique) 
    107      ! 
    108      ! stop if error in size 
    109      IF( imax > SIZE(G0%glamt,1) .OR. jmax > SIZE(G0%glamt,2) .OR. imax <= imin .OR. jmax <= jmin ) THEN                    
    110         WRITE(*,*) 'ERROR ***** bad child grid definition ...' 
    111         WRITE(*,*) 'please check imin,jmin,imax,jmax,jpizoom,jpjzoom values'        
    112         STOP 
    113      ENDIF 
    114      IF( SIZE(G1%nav_lon,1) .NE. nxfin .OR. SIZE(G1%nav_lon,2) .NE. nyfin ) THEN 
    115         WRITE(*,*) 'ERROR ***** bad child coordinates file ...' 
    116         WRITE(*,*) 'please check that child coordinates file has been created with the same namelist' 
    117         STOP 
    118      ENDIF 
    119      ! 
    120      jpj = SIZE(G0%nav_lat,2) 
    121      jpi = SIZE(G0%nav_lat,1)     
    122      !            
    123      ! create logical array masksrc 
    124      ALLOCATE(masksrc(jpi,jpj)) 
    125      WHERE(G0%bathy_level.LE.0)   ;   masksrc = .FALSE.   ; 
    126      ELSEWHERE                    ;   masksrc = .TRUE.    ; 
    127      END WHERE 
    128  
    129      ! change longitudes (from -180:180 to 0:360) 
    130      IF ( Pacifique ) THEN 
    131         WHERE(G0%nav_lon < 0.001)   G0%nav_lon = G0%nav_lon + 360. 
    132      ENDIF 
    133      !              
    134      ALLOCATE(G1%bathy_meter(nxfin,nyfin)) 
    135      ! 
    136      ! compute remapping matrix thanks to SCRIP package 
    137      CALL get_remap_matrix(G0%nav_lat,G1%nav_lat,G0%nav_lon,G1%nav_lon,masksrc,matrix,src_add,dst_add) 
    138      CALL make_remap(G0%bathy_meter,G1%bathy_meter,nxfin,nyfin,matrix,src_add,dst_add)   
    139      DEALLOCATE(masksrc) 
    140      !       
     147     ! 
     148     IF( .NOT. identical_grids ) THEN 
     149        !            
     150        WRITE(*,*) 'Interpolation of high resolution bathymetry on child grid: Bilinear interpolation ...' 
     151        ! 
     152        ! create logical array masksrc 
     153        ALLOCATE(masksrc(SIZE(G0%bathy_meter,1),SIZE(G0%bathy_meter,2))) 
     154        WHERE(G0%bathy_level.LE.0)   ;   masksrc = .FALSE.   ; 
     155        ELSEWHERE                    ;   masksrc = .TRUE.    ; 
     156        END WHERE 
     157        ! 
     158        ! compute remapping matrix thanks to SCRIP package 
     159        CALL get_remap_matrix(G0%nav_lat,G1%nav_lat,G0%nav_lon,G1%nav_lon,masksrc,matrix,src_add,dst_add) 
     160        CALL make_remap(G0%bathy_meter,G1%bathy_meter,nxfin,nyfin,matrix,src_add,dst_add)   
     161        DEALLOCATE(masksrc) 
     162        ! 
     163     ENDIF 
     164      
    141165     ! compute constant bathymetry for Parent-Child bathymetry connection 
    142166     CALL init_constant_bathy(G0%bathy_meter,bathy_fin_constant) 
     
    210234     CALL levels_to_meter(G1) ! needed for domcfg 
    211235     ! 
     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     ! 
    212243     ! store interpolation result in output file 
    213      status = Write_Bathy_level(TRIM(child_level),G1) 
    214      status = write_domcfg(TRIM(child_domcfg),G1) 
     244     IF( TRIM(parent_bathy_level) /= '' )   status = Write_Bathy_level(TRIM(child_level),G1) 
     245     IF( TRIM(parent_bathy_meter) /= '' )   status = Write_Bathy_meter(TRIM(child_meter),G1) 
     246     IF( TRIM(parent_domcfg_out)  /= '' )   status = write_domcfg(TRIM(child_domcfg),G1) 
    215247     !!status = write_domcfg(TRIM(parent_domcfg_out),G0) ! do not activate it 
    216248 
     
    222254     !                                                ! -----------------------------------------------------  
    223255     WRITE(*,*) '*** Second option : new topo or partial steps' 
    224  
    225      ! read fine and coarse grids coordinates file 
    226      status = Read_Coordinates(TRIM(parent_coordinate_file),G0) 
    227      status = Read_Coordinates(TRIM(child_coordinates),G1,Pacifique) 
    228      ! 
    229      ! check error in size 
    230      IF( imax > SIZE(G0%nav_lon,1) .OR. jmax > SIZE(G0%nav_lon,2) .OR. imax <= imin .OR. jmax <= jmin ) THEN                     
    231         WRITE(*,*) 'ERROR ***** bad child grid definition ...' 
    232         WRITE(*,*) 'please check imin,jmin,imax,jmax,jpizoom,jpjzoom values'        
    233         STOP 
    234      ENDIF 
    235      IF( SIZE(G1%nav_lon,1) .NE. nxfin .OR. SIZE(G1%nav_lon,2) .NE. nyfin ) THEN 
    236         WRITE(*,*) 'ERROR ***** bad child coordinates file ...' 
    237         WRITE(*,*) 'please check that child coordinates file has been created with the same namelist' 
    238         STOP 
    239      ENDIF 
    240      ! 
    241      ! === From here on: G0 is the grid associated with the new topography (as gebco or etopo) === 
    242      ! 
    243      ! read bathymetry data set => G0%bathy_meter 
    244      IF( new_topo ) THEN                                       ! read G0%bathy_meter (on a reduced grid) and G1 coordinates 
    245         DEALLOCATE( G0%nav_lon, G0%nav_lat ) 
    246         status = read_bathy_coord(TRIM(elevation_database),G0,G1,Pacifique) 
    247      ELSE                                                      ! read G0%bathy_meter (on the global grid) 
    248         IF( TRIM(parent_bathy_level) /= '' ) THEN 
    249            status = Read_bathy_level(TRIM(parent_bathy_level),G0) 
    250            CALL levels_to_meter(G0) 
    251         ELSEIF( TRIM(parent_bathy_level) == '' .AND. TRIM(parent_bathy_meter) /= '') THEN 
    252            status = read_bathy_meter(TRIM(parent_bathy_meter),G0) 
    253         ENDIF 
    254         IF(Pacifique) THEN 
    255            WHERE(G0%nav_lon < 0.0001)   G0%nav_lon = G0%nav_lon + 360. 
    256         ENDIF 
    257      ENDIF 
    258      !                
     256     ! 
     257     ! === Here: G0 is the grid associated with the new topography (as gebco or etopo) === 
     258     ! 
    259259     ! what type of interpolation for bathymetry 
    260260     IF( type_bathy_interp == 0 ) THEN 
     
    269269     ENDIF 
    270270     ! 
    271      ! 1st allocation of child grid bathy 
    272      ALLOCATE(G1%bathy_meter(nxfin,nyfin)) 
    273      G1%bathy_meter(:,:)=0.                        
    274271     ! 
    275272     ! --------------------------------------------------------------------------------- 
     
    278275     ! ==> It gives G1%bathy_meter from G0%bathy_meter 
    279276     ! --------------------------------------------------------------------------------- 
    280      identical_grids = .FALSE. 
    281  
    282      IF(  SIZE(G0%nav_lat,1) == SIZE(G1%nav_lat,1) .AND. SIZE(G0%nav_lat,2) == SIZE(G1%nav_lat,2) .AND. & 
    283         & SIZE(G0%nav_lon,1) == SIZE(G1%nav_lon,1) .AND. SIZE(G0%nav_lon,2) == SIZE(G1%nav_lon,2) ) THEN 
    284         IF(  MAXVAL( ABS(G0%nav_lat(:,:)- G1%nav_lat(:,:)) ) < 0.0001 .AND. & 
    285            & MAXVAL( ABS(G0%nav_lon(:,:)- G1%nav_lon(:,:)) ) < 0.0001 ) THEN 
    286            WRITE(*,*) '' 
    287            WRITE(*,*) 'same grid between parent and child domains => NO INTERPOLATION' 
    288            WRITE(*,*) '' 
    289            G1%bathy_meter = G0%bathy_meter 
    290            identical_grids = .TRUE. 
    291         ENDIF 
    292      ENDIF 
    293277 
    294278     IF( .NOT. identical_grids ) THEN 
     
    392376           ENDIF 
    393377 
     378           DEALLOCATE(trouble_points) 
     379 
    394380        ENDIF 
    395381        !                                                       ! -----------------------------  
     
    400386           ALLOCATE(bathy_test(nxfin,nyfin)) 
    401387           ! 
    402            !Where(G0%bathy_meter.le.0.00001)  
    403            !  masksrc = .false. 
    404            !ElseWhere 
    405            masksrc = .TRUE. 
    406            !End where                        
     388           WHERE(G0%bathy_level.LE.0)   ;   masksrc = .FALSE.   ; 
     389           ELSEWHERE                    ;   masksrc = .TRUE.    ; 
     390           END WHERE 
    407391           !             
    408392           ! compute remapping matrix thanks to SCRIP package             
     
    446430     ! allocate temporary arrays                   
    447431     IF (.NOT.ASSOCIATED(G0%gdepw_ps))       ALLOCATE(G0%gdepw_ps    (SIZE(G0%bathy_meter,1),SIZE(G0%bathy_meter,2))) 
    448      IF (.NOT.ASSOCIATED(G1%gdepw_ps))       ALLOCATE(G1%gdepw_ps    (SIZE(G1%bathy_meter,1),SIZE(G1%bathy_meter,2)))                   
     432     IF (.NOT.ASSOCIATED(G1%gdepw_ps))       ALLOCATE(G1%gdepw_ps    (SIZE(G1%bathy_meter,1),SIZE(G1%bathy_meter,2))) 
    449433     IF (.NOT.ASSOCIATED(gdepw_ps_interp))   ALLOCATE(gdepw_ps_interp(SIZE(G1%bathy_meter,1),SIZE(G1%bathy_meter,2))) 
    450434     !                        
     
    619603     ! 
    620604     ! update parent grid 
    621      IF(bathy_update)   CALL Update_Parent_Bathy( G0,G1 )  
    622       
    623      IF(bathy_update)   status = Write_Bathy_meter(TRIM(parent_bathy_meter_updated),G0) 
    624      IF(bathy_update)   status = write_domcfg(TRIM(parent_domcfg_updated),G0) 
    625       
     605     IF(bathy_update) THEN 
     606        CALL Update_Parent_Bathy( G0,G1 )  
     607        status = Write_Bathy_meter(TRIM(parent_bathy_meter_updated),G0) 
     608        status = write_domcfg(TRIM(parent_domcfg_updated),G0) 
     609     ENDIF 
     610     ! 
    626611     ! store interpolation result in output file 
    627612     IF( TRIM(parent_bathy_level) /= '' )   status = Write_Bathy_level(TRIM(child_level),G1) 
Note: See TracChangeset for help on using the changeset viewer.