Changeset 12264


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

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

Location:
utils/tools/NESTING/src
Files:
5 edited

Legend:

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

    r12253 r12264  
    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 = 4 ??? 
    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/NESTING/src/agrif_create_bathy.f90

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

    r10025 r12264  
    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 = 4 ??? 
     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/NESTING/src/agrif_readwrite.f90

    r12253 r12264  
    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      ELSE 
    551         ln_zps = 0 
    552         ln_zco = 1 
    553      ENDIF 
     550!!$     ELSE 
     551!!$        ln_zps = 0 
     552!!$        ln_zco = 1 
     553!!$     ENDIF 
    554554 
    555555     CALL zgr_zps( nx, ny, gdept_1d, gdepw_1d, e3t_1d, e3w_1d, zbathy, &        ! input 
     
    866866      zmax = gdepw_1d(N) + e3t_1d(N)        ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_1d(N-1) ) 
    867867      zbathy(:,:) = MIN( zmax ,  zbathy(:,:) )    ! bounded value of bathy (min already set at the end of zgr_bat) 
    868       WHERE( zbathy(:,:) == 0. )   ;   mbathy(:,:) = 0       ! land  : set mbathy to 0 
     868      WHERE( zbathy(:,:) == 0. )   ;   mbathy(:,:) = 0     ! land  : set mbathy to 0 
    869869      ELSEWHERE                    ;   mbathy(:,:) = N-1   ! ocean : initialize mbathy to the max ocean level 
    870870      END WHERE 
    871871 
    872       IF( partial_steps ) THEN 
     872!!$      IF( partial_steps ) THEN 
    873873         ! Compute mbathy for ocean points (i.e. the number of ocean levels) 
    874874         ! find the number of ocean levels such that the last level thickness 
     
    879879            WHERE( 0. < zbathy(:,:) .AND. zbathy(:,:) <= zdepth )   mbathy(:,:) = jk-1 
    880880         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 
     881!!$      ELSE 
     882!!$         DO jk = 1, N 
     883!!$            WHERE( 0. < zbathy(:,:) .AND. zbathy(:,:) >= gdepw_1d(jk) )   mbathy(:,:) = jk-1 
     884!!$         END DO         
     885!!$      ENDIF 
    886886       
    887887      ! Scale factors and depth at T- and W-points 
     
    897897      !!      IF ( ln_isfcav == 1 ) CALL zgr_isf 
    898898      ! 
    899       IF( partial_steps ) THEN 
    900          ! Scale factors and depth at T- and W-points 
    901          ! IF ( ln_isfcav == 0 ) THEN 
    902          DO jj = 1, ny 
    903             DO ji = 1, nx 
    904                ik = mbathy(ji,jj) 
    905                IF( ik > 0 ) THEN               ! ocean point only 
    906                   ! max ocean level case 
    907                   IF( ik == N-1 ) THEN 
    908                      zdepwp = zbathy(ji,jj) 
    909                      ze3tp  = zbathy(ji,jj) - gdepw_1d(ik) 
    910                      ze3wp = 0.5 * e3w_1d(ik) * ( 1. + ( ze3tp/e3t_1d(ik) ) ) 
    911                      e3t_0(ji,jj,ik  ) = ze3tp 
    912                      e3t_0(ji,jj,ik+1) = ze3tp 
    913                      e3w_0(ji,jj,ik  ) = ze3wp 
    914                      e3w_0(ji,jj,ik+1) = ze3tp 
    915                      gdepw_0(ji,jj,ik+1) = zdepwp 
    916                      gdept_0(ji,jj,ik  ) = gdept_1d(ik-1) + ze3wp 
    917                      gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp 
    918                      ! 
    919                   ELSE                         ! standard case 
    920                      IF( zbathy(ji,jj) <= gdepw_1d(ik+1) ) THEN  ;   gdepw_0(ji,jj,ik+1) = zbathy(ji,jj) 
    921                      ELSE                                        ;   gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 
    922                      ENDIF 
    923                      !gm Bug?  check the gdepw_1d 
    924                      !       ... on ik 
    925                      gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) )   & 
    926                         &                             * ((gdept_1d(     ik  ) - gdepw_1d(ik) )   & 
    927                         &                             / ( gdepw_1d(     ik+1) - gdepw_1d(ik) )) 
    928                      e3t_0  (ji,jj,ik) = e3t_1d  (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) )   &  
    929                         &                             / ( gdepw_1d(      ik+1) - gdepw_1d(ik) )  
    930                      e3w_0(ji,jj,ik) = 0.5 * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2. * gdepw_1d(ik) )   & 
    931                         &                     * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 
    932                      !       ... on ik+1 
    933                      e3w_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik) 
    934                      e3t_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik) 
    935                      gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik) 
     899      ! Scale factors and depth at T- and W-points 
     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) 
    936921                  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) 
    937935               ENDIF 
    938             END DO 
     936            ENDIF 
    939937         END DO 
    940          ! 
    941          ! DO jj = 1, ny 
    942          !    DO ji = 1, nx 
    943          !       ik = mbathy(ji,jj) 
    944          !       IF( ik > 0 ) THEN               ! ocean point only 
    945          !          e3tp (ji,jj) = e3t_0(ji,jj,ik) 
    946          !          e3wp (ji,jj) = e3w_0(ji,jj,ik) 
    947          !       ENDIF 
    948          !    END DO 
    949          ! END DO 
    950          ! END IF 
    951       ENDIF 
     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 
    952950      ! 
    953951      ! Scale factors and depth at U-, V-, UW and VW-points 
  • utils/tools/NESTING/src/agrif_types.f90

    r11206 r12264  
    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       ! 
Note: See TracChangeset for help on using the changeset viewer.