Changeset 12412


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

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

Location:
utils/tools_MERGE_2019
Files:
8 edited

Legend:

Unmodified
Added
Removed
  • utils/tools_MERGE_2019/NESTING/src/agrif_connect_topo.f90

    r10025 r12412  
    251251    ! 
    252252    ! 
    253     DO j = 1,jpj 
    254        DO i = 1,jpi 
     253    DO jj = 1,jpj 
     254       DO ji = 1,jpi 
    255255          ! 
    256           IF (Grid%bathy_meter(i,j) .EQ. 0.0 ) THEN 
    257              Grid%bathy_level(i,j)=0. 
     256          IF (Grid%bathy_meter(ji,jj) .EQ. 0.0 ) THEN 
     257             Grid%bathy_level(ji,jj)=0 
    258258          ELSE 
    259259             !  
    260              k1=4 
     260             k1=2  ! clem: minimum levels = 2 ??? 
    261261             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))) EXIT 
     262                IF ((Grid%bathy_meter(ji,jj).GE.gdepw(k1)) & 
     263                     .AND.(Grid%bathy_meter(ji,jj).LE.gdepw(k1+1))) EXIT 
    264264                k1=k1+1 
    265265             END DO 
    266              Grid%bathy_level(i,j)=k1 
     266             Grid%bathy_level(ji,jj)=k1 
    267267             ! 
    268268          ENDIF 
  • 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 
  • utils/tools_MERGE_2019/NESTING/src/agrif_interpolation.f90

    r10381 r12412  
    198198    CHARACTER(*) :: typevar    
    199199 
    200     INTEGER :: nxf,nyf,zx,zy 
    201     INTEGER :: ji,jj,jif,jjf,jic,jjc,jdecx,jdecy 
     200    INTEGER :: nxf,nyf,nxc,nyc,zx,zy 
     201    INTEGER :: ji,jj,jif,jjf,jic,jjc,jic1,jjc1,jdecx,jdecy 
    202202    REAL*8 :: Ax, Bx, Ay, By 
    203203 
    204204    nxf = SIZE(tabout,1) 
    205205    nyf = SIZE(tabout,2) 
     206 
     207    nxc = SIZE(tabin,1) 
     208    nyc = SIZE(tabin,2) 
    206209 
    207210    SELECT CASE(typevar) 
     
    248251          Ay = 1. - By 
    249252 
    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 
    252258       END DO 
    253259    END DO 
  • utils/tools_MERGE_2019/NESTING/src/agrif_partial_steps.f90

    r10025 r12412  
    174174    END DO 
    175175    ! 
     176!!$    IF( partial_steps ) THEN 
    176177    ! Compute bathy_level for ocean points (i.e. the number of ocean levels) 
    177178    ! find the number of ocean levels such that the last level thickness 
    178179    ! is larger than the minimum of e3zps_min and e3zps_rat * e3t (where 
    179180    ! e3t is the reference level thickness    
    180     !      
    181181    DO jk = N-1, 1, -1 
    182182       zdepth = gdepw(jk) + MIN( e3zps_min, e3t(jk)*e3zps_rat ) 
     
    184184          DO ji = 1, jpi 
    185185             IF( 0. < Grid%bathy_meter(ji,jj) .AND. Grid%bathy_meter(ji,jj) <= zdepth ) & 
    186                   Grid%bathy_level(ji,jj) = jk-1 
     186                Grid%bathy_level(ji,jj) = jk-1 
    187187          END DO 
    188188       END DO 
    189189    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 
    191212 
    192213    CALL bathymetry_control(grid%bathy_level) 
     
    675696  SUBROUTINE bathymetry_control(mbathy) 
    676697 
    677     INTEGER ::   i, j, jl            
     698    INTEGER ::   ji, jj, jl            
    678699    INTEGER ::   icompt, ibtest, ikmax           
    679700    REAL*8, DIMENSION(:,:) :: mbathy      
     
    692713    DO jl = 1, 2 
    693714       ! 
    694        DO j = 2, SIZE(mbathy,2)-1 
    695           DO i = 2, SIZE(mbathy,1)-1 
    696  
    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) ) 
    698719             !                
    699              IF( ibtest < mbathy(i,j) ) THEN 
     720             IF( ibtest < mbathy(ji,jj) ) THEN 
    700721                !                   
    701                 WRITE(*,*) 'grid-point(i,j)= ',i,j,'is changed from',mbathy(i,j),' to ', ibtest 
    702                 mbathy(i,j) = ibtest 
     722                WRITE(*,*) 'grid-point(i,j)= ',ji,jj,'is changed from',mbathy(ji,jj),' to ', ibtest 
     723                mbathy(ji,jj) = ibtest 
    703724                icompt = icompt + 1 
    704725                ! 
     
    720741 
    721742    ikmax = 0 
    722     DO j = 1, SIZE(mbathy,2) 
     743    DO jj = 1, SIZE(mbathy,2) 
    723744       DO ji = 1, SIZE(mbathy,1) 
    724           ikmax = MAX( ikmax, NINT(mbathy(i,j)) ) 
     745          ikmax = MAX( ikmax, NINT(mbathy(ji,jj)) ) 
    725746       END DO 
    726747    END DO 
  • utils/tools_MERGE_2019/NESTING/src/agrif_readwrite.f90

    r11231 r12412  
    241241    CALL write_ncdf_var('nav_lon'        ,dimnames,name,Grid%nav_lon    ,'float') 
    242242    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') 
    244244    ! 
    245245    CALL copy_ncdf_att('nav_lon'        ,TRIM(parent_bathy_level),name,MINVAL(Grid%nav_lon),MAXVAL(Grid%nav_lon)) 
     
    545545     ln_sco = 0 
    546546     ln_isfcav = 0 
    547      IF( partial_steps ) THEN 
     547!!$     IF( partial_steps ) THEN 
    548548        ln_zps = 1 
    549549        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 
    556557 
    557558     ! top/bottom levels 
     
    865866      zmax = gdepw_1d(N) + e3t_1d(N)        ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_1d(N-1) ) 
    866867      zbathy(:,:) = MIN( zmax ,  zbathy(:,:) )    ! bounded value of bathy (min already set at the end of zgr_bat) 
    867       WHERE( zbathy(:,:) == 0. )   ;   mbathy(:,:) = 0       ! land  : set mbathy to 0 
     868      WHERE( zbathy(:,:) == 0. )   ;   mbathy(:,:) = 0     ! land  : set mbathy to 0 
    868869      ELSEWHERE                    ;   mbathy(:,:) = N-1   ! ocean : initialize mbathy to the max ocean level 
    869870      END WHERE 
    870871 
    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       
    880887      ! Scale factors and depth at T- and W-points 
    881888      DO jk = 1, N                        ! intitialization to the reference z-coordinate 
     
    891898      ! 
    892899      ! 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) 
    928921                  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) 
    929935               ENDIF 
    930             END DO 
     936            ENDIF 
    931937         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 
    943950      ! 
    944951      ! Scale factors and depth at U-, V-, UW and VW-points 
  • utils/tools_MERGE_2019/NESTING/src/agrif_types.f90

    r11206 r12412  
    192192       ENDIF 
    193193       ! 
     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       ! 
    194199    ELSE 
    195200       ! 
  • utils/tools_MERGE_2019/SIREN/src/create_boundary.F90

    r12080 r12412  
    937937   ELSE 
    938938 
    939       WRITE(cl_errormsg,*) " ERROR : can't find "//TRIM(cl_namelist) 
     939      WRITE(cl_errormsg,*) " ERROR : can't find "//TRIM(cd_namelist) 
    940940      CALL fct_help(cp_myname,cl_errormsg)  
    941941      CALL EXIT(1) 
  • utils/tools_MERGE_2019/SIREN/src/docsrc/5_changeLog.md

    r12080 r12412  
    55# Release 2019-12-03 {#rev_2019-12-03} 
    66## New features 
    7 - M Siren/src/iom_cdf.f90 : 
     7- M src/iom_cdf.f90 : 
    88   - write netcdf file as netcdf4 
    99 
    1010# Release 2019-11-05 {#rev_2019-11-05} 
    1111## New features 
    12 - M Siren/src/function.f90 
     12- M src/function.f90 
    1313- M src/create_bathy.f90 : 
    1414   - add help and version optional arguments 
Note: See TracChangeset for help on using the changeset viewer.