Changeset 866 for palm/trunk/SOURCE/global_min_max.f90

Timestamp:

Mar 28, 2012 6:44:41 AM (12 years ago)

Author:

raasch

Message:

bugfix for timestep calculation in case of Galilei transformation

File:

• palm/trunk/SOURCE/global_min_max.f90 (modified) (6 diffs)

palm/trunk/SOURCE/global_min_max.f90

@@ -1 @@
- SUBROUTINE global_min_max( i1, i2, j1, j2, k1, k2, ar, mode, value, &
+ SUBROUTINE global_min_max( i1, i2, j1, j2, k1, k2, ar, mode, offset, value, &
                             value_ijk, value1, value1_ijk )
 
 !------------------------------------------------------------------------------!
 ! Current revisions:
-! -----------------
+! ------------------
+! new mode "absoff" accounts for an offset in the respective array
 !
 ! Former revisions:
@@ -43 +44 @@
                           fmin_ijk_l(3), value_ijk(3)
     INTEGER, OPTIONAL ::  value1_ijk(3)
-    REAL              ::  value, &
+    REAL              ::  offset, value, &
                           ar(i1:i2,j1:j2,k1:k2)
 #if defined( __ibm )
@@ -62 +63 @@
 !--    Determine the local minimum
        fmin_ijk_l = MINLOC( ar )
-       fmin_ijk_l(1) = i1 + fmin_ijk_l(1) - nbgp    ! MINLOC assumes lowerbound = 1
+       fmin_ijk_l(1) = i1 + fmin_ijk_l(1) - nbgp ! MINLOC assumes lowerbound = 1
        fmin_ijk_l(2) = j1 + fmin_ijk_l(2) - nbgp
        fmin_ijk_l(3) = k1 + fmin_ijk_l(3) - 1
@@ -103 +104 @@
 !--    Determine the local maximum
        fmax_ijk_l = MAXLOC( ar )
-       fmax_ijk_l(1) = i1 + fmax_ijk_l(1) - nbgp    ! MAXLOC assumes lowerbound = 1
+       fmax_ijk_l(1) = i1 + fmax_ijk_l(1) - nbgp ! MAXLOC assumes lowerbound = 1
        fmax_ijk_l(2) = j1 + fmax_ijk_l(2) - nbgp
        fmax_ijk_l(3) = k1 + fmax_ijk_l(3) - 1
@@ -198 +199 @@
 
 !
+!-- Determine absolute maximum of ( array - offset )
+    IF ( mode == 'absoff' )  THEN
+
+!
+!--    Determine the local absolut maximum
+       fmax_l(1)     = 0.0
+       fmax_ijk_l(1) =  i1
+       fmax_ijk_l(2) =  j1
+       fmax_ijk_l(3) =  k1
+       DO  k = k1, k2
+          DO  j = j1, j2
+!
+!--          Attention: the lowest gridpoint is excluded here, because there
+!--          ---------  is no advection at nzb=0 and mode 'absoff' is only
+!--                     used for calculating u,v extrema for CFL-criteria
+             DO  i = i1+1, i2
+                IF ( ABS( ar(i,j,k) - offset ) > fmax_l(1) )  THEN
+                   fmax_l(1) = ABS( ar(i,j,k) - offset )
+                   fmax_ijk_l(1) = i
+                   fmax_ijk_l(2) = j
+                   fmax_ijk_l(3) = k
+                ENDIF
+             ENDDO
+          ENDDO
+       ENDDO
+
+!
+!--    Set a flag in case that the determined value is negative.
+!--    A constant offset has to be subtracted in order to handle the special
+!--    case i=0 correctly
+       IF ( ar(fmax_ijk_l(1),fmax_ijk_l(2),fmax_ijk_l(3)) < 0.0 )  THEN
+          fmax_ijk_l(1) = -fmax_ijk_l(1) - 10
+       ENDIF
+
+#if defined( __parallel )
+       fmax_l(2)  = myid
+       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
+       CALL MPI_ALLREDUCE( fmax_l, fmax, 1, MPI_2REAL, MPI_MAXLOC, comm2d, &
+                           ierr )
+
+!
+!--    Determine the global absolut maximum. Result stored on PE0.
+       id_fmax = fmax(2)
+       IF ( id_fmax /= 0 )  THEN
+          IF ( myid == 0 )  THEN
+             CALL MPI_RECV( fmax_ijk, 3, MPI_INTEGER, id_fmax, 0, comm2d, &
+                            status, ierr )
+          ELSEIF ( myid == id_fmax )  THEN
+             CALL MPI_SEND( fmax_ijk_l, 3, MPI_INTEGER, 0, 0, comm2d, ierr )
+          ENDIF
+       ELSE
+          fmax_ijk = fmax_ijk_l
+       ENDIF
+!
+!--    Send the indices of the just determined absolut maximum to other PEs
+       CALL MPI_BCAST( fmax_ijk, 3, MPI_INTEGER, 0, comm2d, ierr )
+#else
+       fmax(1)  = fmax_l(1)
+       fmax_ijk = fmax_ijk_l
+#endif
+
+    ENDIF
+
+!
 !-- Determine output parameters
     SELECT CASE( mode )
@@ -218 +283 @@
           value1_ijk = fmax_ijk
 
-       CASE( 'abs' )
+       CASE( 'abs', 'absoff' )
 
           value     = fmax(1)