inversion.f90

Go to the documentation of this file.

!********************************************************************
!
!********************************************************************

MODULE inversion

  USE plume_module, ONLY: w0 , r0
  USE variables, ONLY: write_flag

  USE mixture_module, ONLY : mass_flow_rate
  USE rise, ONLY: plume_height, column_regime

  USE inpout, ONLY: write_inversion
  USE rise, ONLY: plumerise

  USE variables, ONLY : wp

  
  IMPLICIT NONE

  REAL(wp) :: opt_value

  REAL(wp) :: opt_height

  REAL(wp) :: opt_mfr

  INTEGER :: opt_regime


  REAL(wp) :: check_height

  SAVE
  
CONTAINS

  !******************************************************************************
  !
  !******************************************************************************
  
  SUBROUTINE invert_height

    USE variables, ONLY: isset

    IMPLICIT NONE

    REAL(wp) :: r_opt , w_opt
    LOGICAL :: search_flag

    IF ( .NOT.isset(w0) ) THEN

       IF ( .NOT.isset(r0) ) THEN

          WRITE(*,*) 'Inversion: Searching for velocity/radius'
          w0 = 100.0_wp
          CALL velocity_radius_search
          
       ELSE

          WRITE(*,*) 'Inversion: Searching for velocity'
          w0 = 100.0_wp
          
          CALL velocity_search(w_opt,search_flag)
          WRITE(*,*) 'Plume_height [m] =',opt_height,search_flag
          WRITE(*,*) 'Velocity [m/s] =',w_opt

          CALL write_inversion(r0,w_opt,opt_mfr,opt_height,search_flag, &
               opt_regime)

         
          write_flag = .true.
          w0 = w_opt
          CALL plumerise
             
       END IF

    ELSE

       IF ( .NOT.isset(r0) ) THEN
          
          WRITE(*,*) 'Inversion: Searching for radius'
          r0 = 50.0_wp
          
          CALL radius_search(r_opt,search_flag)
          WRITE(*,*) 'Best height [m] =',opt_height,search_flag
          WRITE(*,*) 'Radius [m] =',r_opt

          CALL write_inversion(r0,w_opt,opt_mfr,opt_height,search_flag, &
               opt_regime)

          write_flag = .true.
          r0 = r_opt
          CALL plumerise

       ELSE

          WRITE(*,*) 'No Inversion: radius and velocity fixed in input file'
          
          write_flag = .true.
          CALL plumerise

       END IF
          
    END IF
    
  END SUBROUTINE invert_height

  !******************************************************************************
  !
  !******************************************************************************

  SUBROUTINE velocity_radius_search

    USE variables, ONLY: r_min,r_max,n_values
    
    IMPLICIT NONE

    REAL(wp) :: w_opt
    INTEGER :: i
    LOGICAL :: search_flag
        
    WRITE(*,97)
97  FORMAT(1x,'      radius (m) ',1x,' velocity (m/s) ',1x,                     &
         'MER (kg/s)     ',  1x,'plume height (m)',1x,                          &
         ' inversion ',1x,'column regime')

    
    DO i=0,n_values-1

       r0 = r_min * (r_max/r_min)**( i / (n_values-1.0_wp) )
       ! WRITE(*,*) 'r0',r0
       CALL velocity_search(w_opt,search_flag)
       WRITE(*,101) r0,w_opt,opt_mfr,opt_height,search_flag, opt_regime

       CALL write_inversion(r0,w_opt,opt_mfr,opt_height,search_flag,            &
            opt_regime)

    END DO

101 FORMAT(2(2x,f15.8),1(1x,es15.8),1(1x,f15.2)4x,l,7x,i4)

    
  END SUBROUTINE velocity_radius_search

  !******************************************************************************
  !
  !******************************************************************************
  
  SUBROUTINE velocity_search(w_opt,search_flag)
    
    USE variables, ONLY: height_obj, w_min, w_max
    USE variables, ONLY: height_nbl
    USE variables, ONLY : nbl_stop , umbrella_flag

    IMPLICIT none

    REAL(wp),INTENT(OUT) :: w_opt
    LOGICAL,INTENT(OUT) :: search_flag
    REAL(wp) :: w0_init
    REAL(wp) :: w0_0 ,w0_2
    REAL(wp) :: plume_height_0 , plume_height_2
    REAL(wp) :: sign_0 , sign_2
    REAL(wp) :: init_sign , mult_fact

    INTEGER :: iter_interval
    
    write_flag = .false.
    search_flag = .true.

    ! Initial velocity value for the search of the best value
    w0_init = sqrt(w_max*w_min)
    w0 = w0_init

    CALL plumerise

    IF ( ( umbrella_flag ) .OR. ( nbl_stop ) ) THEN
       
       check_height = height_nbl
       
    ELSE
       
       check_height = plume_height
       
    END IF
    
    ! WRITE(*,*) 'first solve',w0,plume_height,INT(column_regime)

    w_opt = w0
    opt_value = abs(check_height-height_obj)
    opt_height = check_height
    opt_mfr = mass_flow_rate
    opt_regime = column_regime
  
    
    IF ( ( check_height .GT. height_obj ) ) THEN

       mult_fact = 1.0_wp/((w_max/w_min)**0.125_wp)
       plume_height_2 = check_height
       
    ELSE

       mult_fact = ((w_max/w_min)**0.125_wp)
       plume_height_0 = check_height

    END IF

    init_sign = check_height-height_obj

    ! loop to search for a velocity interval over which the residual change sign
    search_interval:DO iter_interval=1,4
    
       w0 = (mult_fact**iter_interval)*w0_init
       
       CALL plumerise

       IF ( ( umbrella_flag ) .OR. ( nbl_stop ) ) THEN
          
          check_height = height_nbl
          
       ELSE
          
          check_height = plume_height
          
       END IF

       IF ( abs(check_height-height_obj) .LT. opt_value ) THEN

          w_opt = w0
          opt_value = abs(check_height-height_obj)
          opt_height = check_height
          opt_mfr = mass_flow_rate
          opt_regime = column_regime
          
       END IF

       ! WRITE(*,*) 'search_interval',w0,plume_height,INT(column_regime)

       IF ( (check_height-height_obj)*init_sign .LT. 0.0_wp ) EXIT search_interval
       
    END DO search_interval

    IF ( iter_interval .EQ. 5 ) THEN

       ! WRITE(*,*) 'optimal velocity not found in the interval',w0_init,w0
       w0 = w0_init
       search_flag = .false.
       return

    END IF
    
    init_sign = check_height-height_obj

    IF ( mult_fact .GT. 1.0_wp ) THEN 
    
       w0_2 = w0
       plume_height_2 = check_height
       w0_0 = w0 / mult_fact

    ELSE

       w0_0 = w0
       plume_height_0 = check_height
       w0_2 = w0 / mult_fact

    END IF

    sign_0 = plume_height_0-height_obj
    sign_2 = plume_height_2-height_obj

    search_zero:DO

       w0 = 0.5_wp * ( w0_0 + w0_2 )

       ! WRITE(*,*) 'search_zero',r0,w0,w0_0,w0_2
       
       CALL plumerise

       IF ( ( umbrella_flag ) .OR. ( nbl_stop ) ) THEN
          
          check_height = height_nbl
          
       ELSE
          
          check_height = plume_height
          
       END IF

       IF ( abs(check_height-height_obj) .LT. opt_value ) THEN

          w_opt = w0
          opt_value = abs(check_height-height_obj)
          opt_height = check_height
          opt_mfr = mass_flow_rate
          opt_regime = column_regime

       END IF
       
       !WRITE(*,*) 'plume_height,regime',plume_height,INT(column_regime)
       !WRITE(*,*) 'w0_0,w0_2',w0_0,w0_2
       !WRITE(*,*) 'plume_0,plume_2',plume_height_0,plume_height_2
       !READ(*,*)

       IF ( abs(plume_height_0-plume_height_2) .LT. 1.e-3_wp ) EXIT search_zero
       IF ( abs(check_height-height_obj)/height_obj .LT. 1.e-5_wp ) EXIT search_zero
       IF ( abs(w0_2-w0_0) .LT. 1.e-6_wp ) THEN

          search_flag = .false.
          EXIT search_zero 

       END IF
          
       IF ( (check_height-height_obj)*sign_2 .LT. 0.0_wp ) THEN

          w0_0 = w0
          plume_height_0 = check_height
          sign_0 = plume_height_0-height_obj

      ELSE

          w0_2 = w0
          plume_height_2 = check_height
          sign_2 = plume_height_2-height_obj
          
       END IF

       init_sign = check_height-height_obj
       
    END DO search_zero

    w0 = w0_init
    
  END SUBROUTINE velocity_search

  !******************************************************************************
  !
  !******************************************************************************

  SUBROUTINE radius_search(r_opt,search_flag)
    
    USE variables, ONLY: height_obj , height_nbl
    USE variables, ONLY : nbl_stop , umbrella_flag
    
    IMPLICIT none

    REAL(wp),INTENT(OUT) :: r_opt
    LOGICAL,INTENT(OUT) :: search_flag
    REAL(wp) :: r0_init
    REAL(wp) :: r0_0 ,r0_2
    REAL(wp) :: plume_height_0 , plume_height_2
    REAL(wp) :: sign_0 , sign_2
    REAL(wp) :: init_sign , mult_fact

    INTEGER :: iter_interval
    
    write_flag = .false.
    search_flag = .true.
    
    r0_init = r0

    CALL plumerise
    !WRITE(*,*) 'first solve',r0,plume_height,INT(column_regime)

    r_opt = r0
    
    IF ( ( umbrella_flag ) .OR. ( nbl_stop ) ) THEN

       check_height = height_nbl

    ELSE

       check_height = plume_height

    END IF
    
    opt_value = abs(check_height-height_obj)
    opt_height = check_height

    opt_mfr = mass_flow_rate
    opt_regime = column_regime

    
    IF ( ( check_height .GT. height_obj ) ) THEN

       mult_fact = 0.33_wp
       plume_height_2 = check_height
       
    ELSE

       mult_fact = 3.33_wp
       plume_height_0 = check_height

    END IF

    init_sign = check_height-height_obj

    search_interval:DO iter_interval=1,5
    
       r0 = mult_fact*r0
       
       CALL plumerise

       IF ( ( umbrella_flag ) .OR. ( nbl_stop ) ) THEN
          
          check_height = height_nbl
          
       ELSE
          
          check_height = plume_height
          
       END IF
       
       !WRITE(*,*) 'search interval',r0,plume_height,INT(column_regime)

       IF ( abs(check_height-height_obj) .LT. opt_value ) THEN

          r_opt = r0
          opt_value = abs(check_height-height_obj)
          opt_height = check_height
          opt_mfr = mass_flow_rate
          opt_regime = column_regime

       END IF

       IF ( (check_height-height_obj)*init_sign .LT. 0.0_wp ) EXIT search_interval
       
    END DO search_interval

    !WRITE(*,*) 'iter_interval',iter_interval
    
    IF ( iter_interval .EQ. 6 ) THEN

       !WRITE(*,*) 'optimal velocity not found in the interval',r0_init,r0
       r0 = r0_init
       search_flag = .false.
       RETURN

    END IF
    
    init_sign = plume_height-height_obj

    IF ( mult_fact .GT. 1.0_wp ) THEN 
    
       r0_2 = r0
       plume_height_2 = plume_height
       r0_0 = r0_init

    ELSE

       r0_0 = r0
       plume_height_0 = plume_height
       r0_2 = r0_init


    END IF

    sign_0 = plume_height_0-height_obj
    sign_2 = plume_height_2-height_obj

    search_zero:DO

       r0 = 0.5_wp * ( r0_0 + r0_2 )

       CALL plumerise

       IF ( ( umbrella_flag ) .OR. ( nbl_stop ) ) THEN
          
          check_height = height_nbl
          
       ELSE
          
          check_height = plume_height
          
       END IF
       
       IF ( abs(check_height-height_obj) .LT. opt_value ) THEN

          r_opt = r0
          opt_value = abs(check_height-height_obj)
          opt_height = check_height
          opt_mfr = mass_flow_rate
          opt_regime = column_regime

       END IF
       
       !WRITE(*,*) 'search_zero',r0,check_height,INT(column_regime)
       !WRITE(*,*) 'r0_0,r0_2',r0_0,r0_2
       !WRITE(*,*) 'plume_0,height_obj,plume_2',plume_height_0,height_obj,plume_height_2
       !READ(*,*)

       IF ( abs(plume_height_0-plume_height_2) .LT. 1.e-3_wp ) EXIT search_zero
       IF ( abs(check_height-height_obj)/height_obj .LT. 1.e-5_wp ) EXIT search_zero
       IF ( abs(r0_2-r0_0) .LT. 1.e-6_wp ) THEN

          search_flag = .false.
          EXIT search_zero 

       END IF
          
       IF ( (check_height-height_obj)*sign_2 .LT. 0.0_wp ) THEN

          r0_0 = r0
          plume_height_0 = check_height
          sign_0 = plume_height_0-height_obj

      ELSE

          r0_2 = r0
          plume_height_2 = check_height
          sign_2 = plume_height_2-height_obj
          
       END IF

       init_sign = check_height-height_obj
       
    END DO search_zero

    r0 = r0_init
    
  END SUBROUTINE radius_search
  

  
END MODULE inversion