Changeset 10027

Timestamp:

2018-08-02T17:08:17+02:00 (6 years ago)

Author:

clem

Message:

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

Location:

utils/tools/NESTING/src

Files:

• agrif_create_bathy.f90 (modified) (7 diffs)
• agrif_modutil.f90 (modified) (1 diff)

utils/tools/NESTING/src/agrif_create_bathy.f90

@@ -40 +40 @@
   ! variables declaration
   !      
-  CHARACTER(len=80) :: namelistname
-  CHARACTER*100     :: Childmeter_file,Childlevel_file,Child_coordinates,child_ps, Child_domcfg    
-  LOGICAL :: identical_grids     
-  INTEGER :: i,j,num_links,nb,nbadd,status,narg,iargc     
-  INTEGER :: jpj,jpi,pos,pos2
-  LOGICAL,DIMENSION(:,:),POINTER     :: masksrc => NULL()  
-  INTEGER,DIMENSION(:,:),ALLOCATABLE :: mask_oce,trouble_points
-  INTEGER,DIMENSION(:)  ,POINTER     :: src_add,dst_add => NULL() 
-  REAL*8, DIMENSION(:,:),POINTER     :: matrix,interpdata => NULL()     
-  REAL*8, DIMENSION(:,:),POINTER     :: bathy_fin_constant => NULL()  
-  REAL*8, DIMENSION(:,:),ALLOCATABLE :: bathy_test,vardep,glamhr,gphihr
-  REAL*8, DIMENSION(:)  ,ALLOCATABLE :: vardep1d
-  REAL*8, DIMENSION(:,:),POINTER     :: gdepw_ps_interp => NULL() 
-  REAL*8, DIMENSION(:,:),POINTER     :: save_gdepw,rx,ry,maskedtopo
-  REAL*8  :: Cell_lonmin,Cell_lonmax,Cell_latmin,Cell_latmax,wghts
-  LOGICAL :: Pacifique = .FALSE.
-  INTEGER :: boundary,xpos,ypos,iimin,iimax,jjmax,jjmin
-  INTEGER :: nxhr,nyhr,ji,jj,nbiter
-  INTEGER :: ipt,jpt,iloc,jloc
-  INTEGER, DIMENSION(2) :: i_min,i_max,j_min,j_max
+  CHARACTER(len=80) ::   namelistname
+  CHARACTER*100     ::   Childmeter_file,Childlevel_file,Child_coordinates,child_ps, Child_domcfg    
+  LOGICAL ::   identical_grids     
+  INTEGER ::   nbadd,status,narg,iargc     
+  INTEGER ::   jpj,jpi
+  LOGICAL,DIMENSION(:,:),POINTER     ::   masksrc => NULL()  
+  INTEGER,DIMENSION(:,:),ALLOCATABLE ::   mask_oce,trouble_points
+  INTEGER,DIMENSION(:)  ,POINTER     ::   src_add,dst_add => NULL() 
+  REAL*8, DIMENSION(:,:),POINTER     ::   matrix,interpdata => NULL()     
+  REAL*8, DIMENSION(:,:),POINTER     ::   bathy_fin_constant => NULL()  
+  REAL*8, DIMENSION(:,:),ALLOCATABLE ::   bathy_test,vardep
+  REAL*8, DIMENSION(:)  ,ALLOCATABLE ::   vardep1d
+  REAL*8, DIMENSION(:,:),POINTER     ::   gdepw_ps_interp => NULL() 
+  REAL*8  ::   Cell_lonmin,Cell_lonmax,Cell_latmin,Cell_latmax,wghts
+  LOGICAL ::   Pacifique = .FALSE.
+  INTEGER ::   boundary,iimin,iimax,jjmax,jjmin
+  INTEGER ::   nxhr,nyhr,nxyhr,ji,jj,nbiter
 
   TYPE(Coordinates) :: G0,G1 
@@ -155 +152 @@
      !       
      ! Check errors in bathy
-     DO j=1,nyfin
-        DO i=1,nxfin
-           IF (G1%bathy_level(i,j).LT.0.) THEN
-              PRINT *,'error in ',i,j,G1%bathy_level(i,j)
+     DO jj=1,nyfin
+        DO ji=1,nxfin
+           IF (G1%bathy_level(ji,jj).LT.0.) THEN
+              PRINT *,'error in ',ji,jj,G1%bathy_level(ji,jj)
               STOP
            ENDIF
@@ -178 +175 @@
         DO WHILE (nbadd.NE.0)
            nbadd = 0
-           DO j=2,nyfin-1
-              DO i=2,nxfin-1
-                 IF (G1%bathy_level(i,j).GT.0.) THEN
-                    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
-                       IF (bathy_test(i,j).NE.1.) nbadd = nbadd + 1
-                       bathy_test(i,j)=1.
+           DO jj=2,nyfin-1
+              DO ji=2,nxfin-1
+                 IF (G1%bathy_level(ji,jj).GT.0.) THEN
+                    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
+                       IF (bathy_test(ji,jj).NE.1.) nbadd = nbadd + 1
+                       bathy_test(ji,jj)=1.
                     ENDIF
                  ENDIF
@@ -264 +261 @@
         trouble_points(:,:) = 0
         !                       
-        DO jj = 2,nyfin
-           DO ji = 2,nxfin
+        DO jj = 2, nyfin
+           DO ji = 2, nxfin
               !       
               ! fine grid cell extension               
@@ -325 +322 @@
                  ENDWHERE
                  !
-                 IF( SUM(mask_oce) < 0.5*(nxhr*nyhr) ) THEN ! if more than half of the points are on land then bathy fine = 0
+                 nxyhr = nxhr*nyhr
+                 IF( SUM(mask_oce) < 0.5*(nxyhr) ) THEN ! if more than half of the points are on land then bathy fine = 0
                     G1%bathy_meter(ji,jj) = 0.
                  ELSE
@@ -331 +329 @@
                        G1%bathy_meter(ji,jj) = SUM( vardep(:,:) * mask_oce(:,:) ) / SUM( mask_oce(:,:) )
                     ELSE                                  ! Median average
-                       vardep(:,:) = vardep(:,:) * mask_oce(:,:) - 100000 * ( 1 - mask_oce(:,:) )
-                       ALLOCATE(vardep1d(nxhr*nyhr))
-                       vardep1d = RESHAPE(vardep,(/ nxhr*nyhr /) )
-                       CALL ssort(vardep1d,nxhr*nyhr)
+                       ALLOCATE(vardep1d(nxyhr))
+                       vardep1d = RESHAPE(vardep,(/ nxyhr /) )
+                       !!CALL ssort(vardep1d,nxyhr)
+                       CALL quicksort(vardep1d,1,nxyhr)
                        !
                        ! Calculate median
-                       IF (MOD(SUM(mask_oce),2) .NE. 0) THEN
-                          G1%bathy_meter(ji,jj) = vardep1d( SUM(mask_oce)/2 + 1 )
+                       IF (MOD(nxyhr,2) .NE. 0) THEN
+                          G1%bathy_meter(ji,jj) = vardep1d( nxyhr/2 + 1 )
                        ELSE
-                          G1%bathy_meter(ji,jj) = 0.5 * ( vardep1d(SUM(mask_oce)/2) + vardep1d(SUM(mask_oce)/2+1) )
+                          G1%bathy_meter(ji,jj) = 0.5 * ( vardep1d(nxyhr/2) + vardep1d(nxyhr/2+1) )
                        END IF
                        DEALLOCATE(vardep1d)   
@@ -572 +570 @@
         DO WHILE (nbadd.NE.0)
            nbadd = 0
-           DO j=2,nyfin-1
-              DO i=2,nxfin-1
-                 IF (G1%bathy_meter(i,j).GT.0.) THEN
-                    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
-                       IF (bathy_test(i,j).NE.1.) nbadd = nbadd + 1
-                       bathy_test(i,j)=1.
+           DO jj=2,nyfin-1
+              DO ji=2,nxfin-1
+                 IF (G1%bathy_meter(ji,jj).GT.0.) THEN
+                    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
+                       IF (bathy_test(ji,jj).NE.1.) nbadd = nbadd + 1
+                       bathy_test(ji,jj)=1.
                     ENDIF
 

utils/tools/NESTING/src/agrif_modutil.f90

@@ -49 +49 @@
   END SUBROUTINE ssort
   !
+  !***********************************************************
+  !                    --- quicksort ---
+  ! Author: t-nissie
+  ! License: GPLv3
+  ! Gist: https://gist.github.com/t-nissie/479f0f16966925fa29ea
+  !***********************************************************
+  RECURSIVE SUBROUTINE quicksort(var, first, last)
+     IMPLICIT NONE
+
+     REAL*8, DIMENSION(:), INTENT(inout) :: var
+     INTEGER,              INTENT(in)    :: first, last
+     REAL*8  :: x, t
+     INTEGER :: ji, jj
+     
+     x = var( (first+last) / 2 )
+     ji = first
+     jj = last
+     DO
+        DO WHILE (var(ji) < x)
+           ji=ji+1
+        END DO
+        DO WHILE (x < var(jj))
+           jj=jj-1
+        END DO
+        IF (ji >= jj) EXIT
+        t = var(ji);  var(ji) = var(jj);  var(jj) = t
+        ji=ji+1
+        jj=jj-1
+     END DO
+     IF (first < ji-1) CALL quicksort(var, first, ji-1)
+     IF (jj+1 < last)  CALL quicksort(var, jj+1, last)
+  END SUBROUTINE quicksort
+  
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !************************************************************************