Ignore:
Timestamp:
2018-08-02T17:08:17+02:00 (2 years ago)
Author:
clem
Message:

solve issues with median averages (greatly increase speed and sort out land issues)

File:
1 edited

Legend:

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

    r10025 r10027  
    4040  ! variables declaration 
    4141  !       
    42   CHARACTER(len=80) :: namelistname 
    43   CHARACTER*100     :: Childmeter_file,Childlevel_file,Child_coordinates,child_ps, Child_domcfg     
    44   LOGICAL :: identical_grids      
    45   INTEGER :: i,j,num_links,nb,nbadd,status,narg,iargc      
    46   INTEGER :: jpj,jpi,pos,pos2 
    47   LOGICAL,DIMENSION(:,:),POINTER     :: masksrc => NULL()   
    48   INTEGER,DIMENSION(:,:),ALLOCATABLE :: mask_oce,trouble_points 
    49   INTEGER,DIMENSION(:)  ,POINTER     :: src_add,dst_add => NULL()  
    50   REAL*8, DIMENSION(:,:),POINTER     :: matrix,interpdata => NULL()      
    51   REAL*8, DIMENSION(:,:),POINTER     :: bathy_fin_constant => NULL()   
    52   REAL*8, DIMENSION(:,:),ALLOCATABLE :: bathy_test,vardep,glamhr,gphihr 
    53   REAL*8, DIMENSION(:)  ,ALLOCATABLE :: vardep1d 
    54   REAL*8, DIMENSION(:,:),POINTER     :: gdepw_ps_interp => NULL()  
    55   REAL*8, DIMENSION(:,:),POINTER     :: save_gdepw,rx,ry,maskedtopo 
    56   REAL*8  :: Cell_lonmin,Cell_lonmax,Cell_latmin,Cell_latmax,wghts 
    57   LOGICAL :: Pacifique = .FALSE. 
    58   INTEGER :: boundary,xpos,ypos,iimin,iimax,jjmax,jjmin 
    59   INTEGER :: nxhr,nyhr,ji,jj,nbiter 
    60   INTEGER :: ipt,jpt,iloc,jloc 
    61   INTEGER, DIMENSION(2) :: i_min,i_max,j_min,j_max 
     42  CHARACTER(len=80) ::   namelistname 
     43  CHARACTER*100     ::   Childmeter_file,Childlevel_file,Child_coordinates,child_ps, Child_domcfg     
     44  LOGICAL ::   identical_grids      
     45  INTEGER ::   nbadd,status,narg,iargc      
     46  INTEGER ::   jpj,jpi 
     47  LOGICAL,DIMENSION(:,:),POINTER     ::   masksrc => NULL()   
     48  INTEGER,DIMENSION(:,:),ALLOCATABLE ::   mask_oce,trouble_points 
     49  INTEGER,DIMENSION(:)  ,POINTER     ::   src_add,dst_add => NULL()  
     50  REAL*8, DIMENSION(:,:),POINTER     ::   matrix,interpdata => NULL()      
     51  REAL*8, DIMENSION(:,:),POINTER     ::   bathy_fin_constant => NULL()   
     52  REAL*8, DIMENSION(:,:),ALLOCATABLE ::   bathy_test,vardep 
     53  REAL*8, DIMENSION(:)  ,ALLOCATABLE ::   vardep1d 
     54  REAL*8, DIMENSION(:,:),POINTER     ::   gdepw_ps_interp => NULL()  
     55  REAL*8  ::   Cell_lonmin,Cell_lonmax,Cell_latmin,Cell_latmax,wghts 
     56  LOGICAL ::   Pacifique = .FALSE. 
     57  INTEGER ::   boundary,iimin,iimax,jjmax,jjmin 
     58  INTEGER ::   nxhr,nyhr,nxyhr,ji,jj,nbiter 
    6259 
    6360  TYPE(Coordinates) :: G0,G1  
     
    155152     !        
    156153     ! Check errors in bathy 
    157      DO j=1,nyfin 
    158         DO i=1,nxfin 
    159            IF (G1%bathy_level(i,j).LT.0.) THEN 
    160               PRINT *,'error in ',i,j,G1%bathy_level(i,j) 
     154     DO jj=1,nyfin 
     155        DO ji=1,nxfin 
     156           IF (G1%bathy_level(ji,jj).LT.0.) THEN 
     157              PRINT *,'error in ',ji,jj,G1%bathy_level(ji,jj) 
    161158              STOP 
    162159           ENDIF 
     
    178175        DO WHILE (nbadd.NE.0) 
    179176           nbadd = 0 
    180            DO j=2,nyfin-1 
    181               DO i=2,nxfin-1 
    182                  IF (G1%bathy_level(i,j).GT.0.) THEN 
    183                     IF (MAX(bathy_test(i,j+1),bathy_test(i,j-1),bathy_test(i-1,j),bathy_test(i+1,j)).EQ.1) THEN 
    184                        IF (bathy_test(i,j).NE.1.) nbadd = nbadd + 1 
    185                        bathy_test(i,j)=1. 
     177           DO jj=2,nyfin-1 
     178              DO ji=2,nxfin-1 
     179                 IF (G1%bathy_level(ji,jj).GT.0.) THEN 
     180                    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 
     181                       IF (bathy_test(ji,jj).NE.1.) nbadd = nbadd + 1 
     182                       bathy_test(ji,jj)=1. 
    186183                    ENDIF 
    187184                 ENDIF 
     
    264261        trouble_points(:,:) = 0 
    265262        !                        
    266         DO jj = 2,nyfin 
    267            DO ji = 2,nxfin 
     263        DO jj = 2, nyfin 
     264           DO ji = 2, nxfin 
    268265              !        
    269266              ! fine grid cell extension                
     
    325322                 ENDWHERE 
    326323                 ! 
    327                  IF( SUM(mask_oce) < 0.5*(nxhr*nyhr) ) THEN ! if more than half of the points are on land then bathy fine = 0 
     324                 nxyhr = nxhr*nyhr 
     325                 IF( SUM(mask_oce) < 0.5*(nxyhr) ) THEN ! if more than half of the points are on land then bathy fine = 0 
    328326                    G1%bathy_meter(ji,jj) = 0. 
    329327                 ELSE 
     
    331329                       G1%bathy_meter(ji,jj) = SUM( vardep(:,:) * mask_oce(:,:) ) / SUM( mask_oce(:,:) ) 
    332330                    ELSE                                  ! Median average 
    333                        vardep(:,:) = vardep(:,:) * mask_oce(:,:) - 100000 * ( 1 - mask_oce(:,:) ) 
    334                        ALLOCATE(vardep1d(nxhr*nyhr)) 
    335                        vardep1d = RESHAPE(vardep,(/ nxhr*nyhr /) ) 
    336                        CALL ssort(vardep1d,nxhr*nyhr) 
     331                       ALLOCATE(vardep1d(nxyhr)) 
     332                       vardep1d = RESHAPE(vardep,(/ nxyhr /) ) 
     333                       !!CALL ssort(vardep1d,nxyhr) 
     334                       CALL quicksort(vardep1d,1,nxyhr) 
    337335                       ! 
    338336                       ! Calculate median 
    339                        IF (MOD(SUM(mask_oce),2) .NE. 0) THEN 
    340                           G1%bathy_meter(ji,jj) = vardep1d( SUM(mask_oce)/2 + 1 ) 
     337                       IF (MOD(nxyhr,2) .NE. 0) THEN 
     338                          G1%bathy_meter(ji,jj) = vardep1d( nxyhr/2 + 1 ) 
    341339                       ELSE 
    342                           G1%bathy_meter(ji,jj) = 0.5 * ( vardep1d(SUM(mask_oce)/2) + vardep1d(SUM(mask_oce)/2+1) ) 
     340                          G1%bathy_meter(ji,jj) = 0.5 * ( vardep1d(nxyhr/2) + vardep1d(nxyhr/2+1) ) 
    343341                       END IF 
    344342                       DEALLOCATE(vardep1d)    
     
    572570        DO WHILE (nbadd.NE.0) 
    573571           nbadd = 0 
    574            DO j=2,nyfin-1 
    575               DO i=2,nxfin-1 
    576                  IF (G1%bathy_meter(i,j).GT.0.) THEN 
    577                     IF (MAX(bathy_test(i,j+1),bathy_test(i,j-1),bathy_test(i-1,j),bathy_test(i+1,j)).EQ.1) THEN 
    578                        IF (bathy_test(i,j).NE.1.) nbadd = nbadd + 1 
    579                        bathy_test(i,j)=1. 
     572           DO jj=2,nyfin-1 
     573              DO ji=2,nxfin-1 
     574                 IF (G1%bathy_meter(ji,jj).GT.0.) THEN 
     575                    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 
     576                       IF (bathy_test(ji,jj).NE.1.) nbadd = nbadd + 1 
     577                       bathy_test(ji,jj)=1. 
    580578                    ENDIF 
    581579 
Note: See TracChangeset for help on using the changeset viewer.