text:
!
! CRTM_Atmosphere
!
! Module for adding layers to the CRTM atmosphere structure as required.
!
!
! CREATION HISTORY:
!       Written by:     Paul van Delst, 29-Oct-2007
!                       paul.vandelst@noaa.gov
!

MODULE CRTM_Atmosphere

  ! -----------------
  ! Environment setup
  ! -----------------
  ! Module use
  USE Type_Kinds            , ONLY: fp
  USE Message_Handler       , ONLY: SUCCESS, FAILURE, WARNING, Display_Message
  USE CRTM_Parameters       , ONLY: ZERO, ONE, POINT_5, SET, &
                                    TOA_PRESSURE           , &
                                    MINIMUM_ABSORBER_AMOUNT, &
                                    WATER_CONTENT_THRESHOLD
  USE CRTM_Atmosphere_Define, ONLY: CRTM_Atmosphere_type    , &
                                    OPERATOR(==), &
                                    OPERATOR(+), &
                                    CRTM_Atmosphere_Associated, &
                                    CRTM_Atmosphere_Create, &
                                    CRTM_Atmosphere_AddLayerCopy, &
                                    CRTM_Atmosphere_Zero
  USE CRTM_Model_Profiles   , ONLY: MODEL_LEVEL_PRESSURE, &
                                    CRTM_Get_Model_Profile
  ! ...Internal variable definition module
  USE iAtm_Define,            ONLY: iAtm_type      , &
                                    iAtm_Associated, &
                                    iAtm_Create    , &
                                    iAtm_Destroy

  ! Disable implicit typing
  IMPLICIT NONE

  ! ------------
  ! Visibilities
  ! ------------
  ! Everything private by default
  PRIVATE
  ! Module procedures
  PUBLIC :: CRTM_Atmosphere_Coverage
  PUBLIC :: CRTM_Atmosphere_IsClear
  PUBLIC :: CRTM_Atmosphere_IsFractional
  PUBLIC :: CRTM_Atmosphere_IsOvercast

  PUBLIC :: CRTM_Atmosphere_AddLayers
  PUBLIC :: CRTM_Atmosphere_AddLayers_TL
  PUBLIC :: CRTM_Atmosphere_AddLayers_AD

  PUBLIC :: CRTM_Atmosphere_ClearSkyCopy
  PUBLIC :: CRTM_Atmosphere_ClearSkyCopy_TL
  PUBLIC :: CRTM_Atmosphere_ClearSkyCopy_AD
  ! iAtm entities
  ! ...Structure
  PUBLIC :: iAtm_type
  ! ...Procedures
  PUBLIC :: iAtm_Associated
  PUBLIC :: iAtm_Create
  PUBLIC :: iAtm_Destroy

  ! -----------------
  ! Module parameters
  ! -----------------
  ! Message string length
  INTEGER, PARAMETER :: ML = 256

  ! The cloud coverage type
  INTEGER, PARAMETER :: CLEAR      = -1
  INTEGER, PARAMETER :: FRACTIONAL =  0
  INTEGER, PARAMETER :: OVERCAST   =  1
!  CHARACTER(*), PARAMETER, DIMENSION( -1:1 ) :: &
!    CLOUD_COVERAGE_NAME = [ 'Clear sky          ', &
!                            'Fractional coverage', &
!                            'Overcast           '  ]

CONTAINS

!################################################################################
!################################################################################
!##                                                                            ##
!##                         ## PUBLIC MODULE ROUTINES ##                       ##
!##                                                                            ##
!################################################################################
!################################################################################

  FUNCTION CRTM_Atmosphere_IsClear(coverage_flag) RESULT(is_clear)
    INTEGER, INTENT(IN) :: coverage_flag
    LOGICAL :: is_clear
    is_clear = coverage_flag == CLEAR
  END FUNCTION CRTM_Atmosphere_IsClear

  FUNCTION CRTM_Atmosphere_IsFractional(coverage_flag) RESULT(is_fractional)
    INTEGER, INTENT(IN) :: coverage_flag
    LOGICAL :: is_fractional
    is_fractional = coverage_flag == FRACTIONAL
  END FUNCTION CRTM_Atmosphere_IsFractional

  FUNCTION CRTM_Atmosphere_IsOvercast(coverage_flag) RESULT(is_overcast)
    INTEGER, INTENT(IN) :: coverage_flag
    LOGICAL :: is_overcast
    is_overcast = coverage_flag == OVERCAST
  END FUNCTION CRTM_Atmosphere_IsOvercast

!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!   CRTM_Atmosphere_Coverage
!
! PURPOSE:
!   Function to determine the cloud coverage type for an input
!   atmosphere.
!
! CALLING SEQUENCE:
!   coverage_flag = CRTM_Atmosphere_Coverage( atm )  ! Input
!
! INPUTS:
!   atm:            Atmosphere structure for which the coverage type is
!                   to be determined.
!                   UNITS:      N/A
!                   TYPE:       CRTM_Atmosphere_type
!                   DIMENSION:  Scalar
!                   ATTRIBUTES: INTENT(IN)
!
! FUNCTION RESULT:
!   coverage_flag:  An integer defining the coverage type. Valid
!                   parameterised values are:
!                     CLEAR
!                     FRACTIONAL
!                     OVERCAST
!                   UNITS:      N/A
!                   TYPE:       INTEGER
!                   DIMENSION:  Scalar
!
!:sdoc-:
!--------------------------------------------------------------------------------
   FUNCTION CRTM_Atmosphere_Coverage(atm) RESULT(coverage_flag)
    ! Arguments
    TYPE(CRTM_Atmosphere_type), INTENT(IN) :: atm
    ! Function result
    INTEGER :: coverage_flag
    ! Local parameters
    REAL(fp), PARAMETER :: MIN_COVERAGE_THRESHOLD = 1.0e-06_fp
    REAL(fp), PARAMETER :: MAX_COVERAGE_THRESHOLD = ONE - MIN_COVERAGE_THRESHOLD
    ! Local variables
    LOGICAL :: cloudy_layer_mask(atm%n_Layers)
    INTEGER :: idx(atm%n_Layers)
    INTEGER :: n, nc, k

    ! Default clear
    coverage_flag = CLEAR
    IF ( atm%n_Clouds == 0 ) RETURN

    ! Check each cloud separately
    Cloud_Loop: DO n = 1, atm%n_Clouds

      ! Determine if there are ANY cloudy layers
      cloudy_layer_mask = atm%Cloud(n)%Water_Content > WATER_CONTENT_THRESHOLD
      nc = COUNT(cloudy_layer_mask)
      IF ( nc == 0 ) CYCLE Cloud_Loop

      ! Get the indices of those cloudy layers
      idx(1:nc) = PACK([(k, k=1,atm%Cloud(n)%n_Layers)], cloudy_layer_mask)

      ! Check for ANY fractional coverage
      ! ??? How to do this without the loop ???
      DO k = 1, nc
!       IF ( (atm%Cloud_Fraction(idx(k)) > MIN_COVERAGE_THRESHOLD) .AND. &
!            (atm%Cloud_Fraction(idx(k)) < MAX_COVERAGE_THRESHOLD) ) THEN
        IF ( (atm%Cloud_Fraction(idx(k)) > MIN_COVERAGE_THRESHOLD) ) THEN
          coverage_flag = FRACTIONAL
          RETURN
        END IF
      END DO

!     ! Check for ALL totally clear or totally cloudy
!     IF ( ALL(atm%Cloud_Fraction(idx(1:nc)) < MIN_COVERAGE_THRESHOLD) .OR. &
!          ALL(atm%Cloud_Fraction(idx(1:nc)) > MAX_COVERAGE_THRESHOLD) ) coverage_flag = OVERCAST

    END DO Cloud_Loop

  END FUNCTION CRTM_Atmosphere_Coverage

!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       CRTM_Atmosphere_AddLayers
!
! PURPOSE:
!       Function to copy an atmosphere structure and adding extra layers from
!       climatology as required to supplement the upper atmosphere profile data.
!
! CALLING SEQUENCE:
!       Error_Status = CRTM_Atmosphere_AddLayers( Atm_In , &  ! Input
!                                                 Atm_Out  )  ! Output
!
! INPUTS:
!       Atm_In:          Atmosphere structure that is to be supplemented
!                        or copied.
!                        UNITS:      N/A
!                        TYPE:       CRTM_Atmosphere_type
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(IN)
!
! OUTPUTS:
!       Atm_Out:         Copy of the input atmosphere structure with extra upper
!                        atmosphere layers added as required.
!                        UNITS:      N/A
!                        TYPE:       CRTM_Atmosphere_type
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(IN OUT)
!
! FUNCTION RESULT:
!       Error_Status:    The return value is an integer defining the error status.
!                        The error codes are defined in the Message_Handler module.
!                        If == SUCCESS the operation was successful
!                           == FAILURE an error occurred
!                        UNITS:      N/A
!                        TYPE:       INTEGER
!                        DIMENSION:  Scalar
!
!:sdoc-:
!--------------------------------------------------------------------------------

  FUNCTION CRTM_Atmosphere_AddLayers( &
    Atm_In , &  ! Input
    Atm_Out) &  ! Output
  RESULT( err_stat )
    ! Arguments
    TYPE(CRTM_Atmosphere_type), INTENT(IN)     :: Atm_In
    TYPE(CRTM_Atmosphere_type), INTENT(IN OUT) :: Atm_Out
    ! Function result
    INTEGER :: err_stat
    ! Local parameters
    CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'CRTM_Atmosphere_AddLayers'
    ! Local variables
    CHARACTER(ML) :: msg
    INTEGER :: i, j, k, n
    TYPE(iAtm_type) :: iAtm

    ! Set up
    err_stat = SUCCESS

    ! If extra layers are NOT needed,
    ! then simply copy the structure
    IF ( Atm_In%Level_Pressure(0) <= TOA_PRESSURE) THEN
      Atm_Out = Atm_In
      IF ( .NOT. CRTM_Atmosphere_Associated( Atm_Out ) ) THEN
        err_stat = FAILURE
        msg = 'Error assigning Atmosphere structure with NO extra layers'
        CALL Display_Message( ROUTINE_NAME, msg, err_stat )
      END IF
      RETURN
    END IF

    ! Determine the number of extra layers required
    n = Extra_Layers( Atm_In )
    IF ( n < 1 ) THEN
      err_stat = FAILURE
      msg = 'Error determining extra layer count'
      CALL Display_Message( ROUTINE_NAME, msg, err_stat )
      RETURN
    END IF

    ! Allocate the internal variable structure
    CALL iAtm_Create( iAtm, n, Atm_In%n_Absorbers )
    IF ( .NOT. iAtm_Associated( iAtm ) ) THEN
      err_stat = FAILURE
      msg = 'Error allocating iAtm internal structure'
      CALL Display_Message( ROUTINE_NAME, msg, err_stat )
      RETURN
    END IF

    ! Get the extra layer profiles
    CALL CRTM_Get_Model_Profile( Atm_In%Absorber_Id, &
                                 iAtm%pl, iAtm%tl, iAtm%al, &
                                 Model=Atm_In%Climatology )

    ! First interpolate the extra levels to the user top pressure
    ! replacing the model data at that array index
    CALL Interp_LPoly( Atm_In%Level_Pressure(0), iAtm%pl(n-1:n), iAtm%ilpoly )
    iAtm%plint_save = Atm_In%Level_Pressure(0)
    iAtm%pln_save   = iAtm%pl(n)
    iAtm%pl(n)      = iAtm%plint_save
    CALL Interp_Linear( iAtm%ilpoly, iAtm%tl(n-1:n), iAtm%tlint_save )
    iAtm%tln_save = iAtm%tl(n)
    iAtm%tl(n)    = iAtm%tlint_save
    DO j = 1, Atm_In%n_Absorbers
      CALL Interp_Linear( iAtm%ilpoly, iAtm%al(n-1:n,j), iAtm%alint_save(j) )
      iAtm%aln_save(j) = iAtm%al(n,j)
      iAtm%al(n,j)     = iAtm%alint_save(j)
    END DO

    ! Now compute the model profile layer averages
    DO k = 1, n
      CALL Layer_P(iAtm%pl(k-1:k), iAtm%p(k))
      CALL Layer_X(iAtm%tl(k-1:k), iAtm%t(k))
    END DO
    DO j = 1, Atm_In%n_Absorbers
      DO k = 1, n
        CALL Layer_X(iAtm%al(k-1:k,j), iAtm%a(k,j))
      END DO
    END DO

    ! Now, extrapolate user layer profile to get the "layer 0" value and
    ! use it to shift the model profile to avoid a discontinuity at p(n)
    CALL Interp_LPoly( iAtm%p(n), Atm_In%Pressure(1:2), iAtm%elpoly )
    CALL Shift_Profile( iAtm%elpoly, Atm_In%Temperature(1:2), iAtm%t )
    DO j = 1, Atm_In%n_Absorbers
      CALL Shift_Profile( iAtm%elpoly, Atm_In%Absorber(1:2,j), iAtm%a(:,j) )
    END DO

    ! Make sure the absorber amounts are not negative.
    ! (Is a further, more physical, check needed here?)
    iAtm%a_save = iAtm%a
    WHERE (iAtm%a_save < ZERO) iAtm%a = MINIMUM_ABSORBER_AMOUNT

    ! Copy over the atmosphere structure with extra layers
    atm_out = CRTM_Atmosphere_AddLayerCopy( Atm_In, n )
    IF ( .NOT. CRTM_Atmosphere_Associated( atm_out ) ) THEN
      err_stat = FAILURE
      msg = 'Error copying Atmosphere structure with extra layers'
      CALL Display_Message( ROUTINE_NAME, msg, err_stat )
      RETURN
    END IF

    ! Slot the added layers into the output atmosphere structure
    ! Note: Cloud and Aerosol assignments not really needed (the
    !       zeroing is handled by the structure allocation) since
    !       at TOA, typically, there are not any clouds and/or
    !       aerosols.
    ! ...Profile
    Atm_Out%Level_Pressure(0:n) = iAtm%pl
    Atm_Out%Pressure(1:n)       = iAtm%p
    Atm_Out%Temperature(1:n)    = iAtm%t
    Atm_Out%Cloud_Fraction(1:n) = ZERO
    DO j = 1, Atm_Out%n_Absorbers
      Atm_Out%Absorber(1:n,j)   = iAtm%a(:,j)
    END DO
    ! ...Clouds
    IF ( Atm_In%n_Clouds > 0 ) THEN
      DO i = 1, Atm_In%n_Clouds
        Atm_Out%Cloud(i)%Effective_Radius(1:n)   = ZERO
        Atm_Out%Cloud(i)%Effective_Variance(1:n) = ZERO
        Atm_Out%Cloud(i)%Water_Content(1:n)      = ZERO
      END DO
    END IF
    ! ...Aerosols
    IF ( Atm_In%n_Aerosols > 0 ) THEN
      DO i = 1, Atm_In%n_Aerosols
        Atm_Out%Aerosol(i)%Effective_Radius(1:n) = ZERO
        Atm_Out%Aerosol(i)%Concentration(1:n)    = ZERO
      END DO
    END IF

    ! Clean up
    CALL iAtm_Destroy( iAtm )

  END FUNCTION CRTM_Atmosphere_AddLayers

!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       CRTM_Atmosphere_AddLayers_TL
!
! PURPOSE:
!       Function to copy a tangent-linear atmosphere structure and add extra
!       layers as required to supplement the upper atmosphere profile data.
!
! CALLING SEQUENCE:
!       Error_Status = CRTM_Atmosphere_AddLayers_TL( Atm_In    , &  ! FWD Input
!                                                    Atm_In_TL , &  ! TL  Input
!                                                    Atm_Out_TL  )  ! TL  Output
!
! INPUTS:
!       Atm_In:          Forward model atmosphere structure.
!                        UNITS:      N/A
!                        TYPE:       CRTM_Atmosphere_type
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(IN)
!
!       Atm_In_TL:       Tangent-linear model atmosphere structure that is
!                        to be copied.
!                        UNITS:      N/A
!                        TYPE:       CRTM_Atmosphere_type
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(IN)
!
! OUTPUTS:
!       Atm_Out_TL:      Copy of the input tangent-linear atmosphere structure
!                        with extra upper atmosphere layers added as required.
!                        Note that the tangent-linear values of the added layers
!                        is *always* zero.
!                        UNITS:      N/A
!                        TYPE:       CRTM_Atmosphere_type
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(IN OUT)
!
! FUNCTION RESULT:
!       Error_Status:    The return value is an integer defining the error status.
!                        The error codes are defined in the Message_Handler module.
!                        If == SUCCESS the operation was successful
!                           == FAILURE an error occurred
!                        UNITS:      N/A
!                        TYPE:       INTEGER
!                        DIMENSION:  Scalar
!
!:sdoc-:
!--------------------------------------------------------------------------------

  FUNCTION CRTM_Atmosphere_AddLayers_TL( &
    Atm_In    , &  ! FWD Input
    Atm_In_TL , &  ! TL  Input
    Atm_Out_TL) &  ! TL  Output
  RESULT( err_stat )
    ! Arguments
    TYPE(CRTM_Atmosphere_type), INTENT(IN)     :: Atm_In
    TYPE(CRTM_Atmosphere_type), INTENT(IN)     :: Atm_In_TL
    TYPE(CRTM_Atmosphere_type), INTENT(IN OUT) :: Atm_Out_TL
    ! Function result
    INTEGER :: err_stat
    ! Local parameters
    CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'CRTM_Atmosphere_AddLayers_TL'
    ! Local variables
    CHARACTER(ML) :: msg
    INTEGER :: n

    ! Set up
    err_stat = SUCCESS

    ! If extra layers are NOT needed,
    ! then simply copy the structure
    IF ( Atm_In%Level_Pressure(0) <= TOA_PRESSURE) THEN
      Atm_Out_TL = Atm_In_TL
      IF ( .NOT. CRTM_Atmosphere_Associated( Atm_Out_TL ) ) THEN
        err_stat = FAILURE
        msg = 'Error assigning Atmosphere structure with NO extra layers'
        CALL Display_Message( ROUTINE_NAME, msg, err_stat )
      END IF
      RETURN
    END IF

    ! Determine how many extra layers are needed
    n = Extra_Layers( Atm_In )
    IF ( n < 1 ) THEN
      err_stat = FAILURE
      msg = 'Error determining extra layer count'
      CALL Display_Message( ROUTINE_NAME, msg, err_stat )
      RETURN
    END IF

    ! Copy over the atmosphere structure with extra layers
    ! (which will be zero by definition)
    atm_out_TL = CRTM_Atmosphere_AddLayerCopy( Atm_In_TL, n )
    IF ( .NOT. CRTM_Atmosphere_Associated( atm_out_TL ) ) THEN
      err_stat = FAILURE
      msg = 'Error copying Atmosphere structure with extra layers'
      CALL Display_Message( ROUTINE_NAME, msg, err_stat )
      RETURN
    END IF

  END FUNCTION CRTM_Atmosphere_AddLayers_TL

!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       CRTM_Atmosphere_AddLayers_AD
!
! PURPOSE:
!       Function to copy back an adjoint atmosphere structure removing added
!       extra layers as were required to supplement the upper atmosphere
!       profile data.
!
! CALLING SEQUENCE:
!       Error_Status = CRTM_Atmosphere_AddLayers_AD( Atm_In    , &  ! FWD Input
!                                                    Atm_Out_AD, &  ! AD  Input
!                                                    Atm_In_AD   )  ! AD  Output
!
! INPUTS:
!       Atm_In:          Forward model atmosphere structure.
!                        UNITS:      N/A
!                        TYPE:       CRTM_Atmosphere_type
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(IN)
!
!       Atm_Out_AD:      Adjoint atmosphere structure that contains the added
!                        extra layers.
!                        ** SET TO ZERO ON EXIT **
!                        UNITS:      N/A
!                        TYPE:       CRTM_Atmosphere_type
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(IN OUT)
!
! OUTPUTS:
!       Atm_In_AD:       Adjoint atmosphere structure at the original, user
!                        specified layering.
!                        ** MUST HAVE VALUE ON ENTRY **
!                        UNITS:      N/A
!                        TYPE:       CRTM_Atmosphere_type
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(IN OUT)
!
! FUNCTION RESULT:
!       Error_Status:    The return value is an integer defining the error status.
!                        The error codes are defined in the Message_Handler module.
!                        If == SUCCESS the operation was successful
!                           == FAILURE an error occurred
!                        UNITS:      N/A
!                        TYPE:       INTEGER
!                        DIMENSION:  Scalar
!
!:sdoc-:
!--------------------------------------------------------------------------------

  FUNCTION CRTM_Atmosphere_AddLayers_AD( &
    Atm_In    , &  ! FWD Input
    Atm_Out_AD, &  ! AD  Input
    Atm_In_AD ) &  ! AD  Output
  RESULT( err_stat )
    ! Arguments
    TYPE(CRTM_Atmosphere_type), INTENT(IN)     :: Atm_In
    TYPE(CRTM_Atmosphere_type), INTENT(IN OUT) :: Atm_Out_AD
    TYPE(CRTM_Atmosphere_type), INTENT(IN OUT) :: Atm_In_AD
    ! Function result
    INTEGER :: err_stat
    ! Local parameters
    CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'CRTM_Atmosphere_AddLayers_AD'
    ! Local variables
    CHARACTER(ML) :: msg
    INTEGER :: i, j, n, no, nt

    ! Set up
    err_stat = SUCCESS

    ! If extra layers are NOT needed, then simply perform
    ! the adjoint sum. Remember the TL form is
    !   Atm_Out_TL = Atm_In_TL
    ! so the adjoint form is
    !   Atm_In_AD  = Atm_In_AD + Atm_Out_AD
    !   Atm_Out_AD = ZERO
    IF ( Atm_In%Level_Pressure(0) <= TOA_PRESSURE) THEN
      Atm_In_AD = Atm_In_AD + Atm_Out_AD
      CALL CRTM_Atmosphere_Zero( Atm_Out_AD )
      RETURN
    END IF

    ! Determine how many extra layers have been used
    n = Extra_Layers( Atm_In )
    IF ( n < 1 ) THEN
      err_stat = FAILURE
      msg = 'Error determining extra layer count'
      CALL Display_Message( ROUTINE_NAME, msg, err_stat )
      RETURN
    END IF

    ! Perform the adjoint summations
    ! This is the adjoint equivalent of the TL Assign_Atmosphere
    no = Atm_In_AD%n_Layers
    nt = n + no
    ! ...Aerosols
    IF ( Atm_In_AD%n_Aerosols > 0 ) THEN
      DO i = 1, Atm_In_AD%n_Aerosols
        Atm_In_AD%Aerosol(i)%Concentration(1:no) = Atm_In_AD%Aerosol(i)%Concentration(1:no) + &
                                                   Atm_Out_AD%Aerosol(i)%Concentration(n+1:nt)
        Atm_In_AD%Aerosol(i)%Effective_Radius(1:no) = Atm_In_AD%Aerosol(i)%Effective_Radius(1:no) + &
                                                      Atm_Out_AD%Aerosol(i)%Effective_Radius(n+1:nt)
        Atm_In_AD%Aerosol(i)%Type = Atm_Out_AD%Aerosol(i)%Type
      END DO
    END IF
    ! ...Clouds
    IF ( Atm_In_AD%n_Clouds > 0 ) THEN
      DO i = 1, Atm_In_AD%n_Clouds
        Atm_In_AD%Cloud(i)%Water_Content(1:no) = Atm_In_AD%Cloud(i)%Water_Content(1:no) + &
                                                 Atm_Out_AD%Cloud(i)%Water_Content(n+1:nt)
        Atm_In_AD%Cloud(i)%Effective_Variance(1:no) = Atm_In_AD%Cloud(i)%Effective_Variance(1:no) + &
                                                      Atm_Out_AD%Cloud(i)%Effective_Variance(n+1:nt)
        Atm_In_AD%Cloud(i)%Effective_Radius(1:no) = Atm_In_AD%Cloud(i)%Effective_Radius(1:no) + &
                                                    Atm_Out_AD%Cloud(i)%Effective_Radius(n+1:nt)
        Atm_In_AD%Cloud(i)%Type = Atm_Out_AD%Cloud(i)%Type
      END DO
    END IF
    ! ...Absorber data
    DO j = 1, Atm_In_AD%n_Absorbers
      Atm_In_AD%Absorber(1:no,j) = Atm_In_AD%Absorber(1:no,j) + Atm_Out_AD%Absorber(n+1:nt,j)
    END DO
    ! ...Temperature data
    Atm_In_AD%Temperature(1:no) = Atm_In_AD%Temperature(1:no) + Atm_Out_AD%Temperature(n+1:nt)
    ! ...Pressure data
    Atm_In_AD%Pressure(1:no)       = Atm_In_AD%Pressure(1:no) + Atm_Out_AD%Pressure(n+1:nt)
    Atm_In_AD%Level_Pressure(0:no) = Atm_In_AD%Level_Pressure(0:no) + Atm_Out_AD%Level_Pressure(n:nt)
    ! ...Cloud fraction data
    Atm_In_AD%Cloud_Fraction(1:no) = Atm_In_AD%Cloud_Fraction(1:no) + &
                                     Atm_Out_AD%Cloud_Fraction(n+1:nt)

    ! Zero the output atmosphere structure
    CALL CRTM_Atmosphere_Zero( Atm_Out_AD )

  END FUNCTION CRTM_Atmosphere_AddLayers_AD

!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       CRTM_Atmosphere_ClearSkyCopy
!
! PURPOSE:
!       Function to copy an instance of the CRTM Atmosphere object
!       but without the clouds included.
!
! CALLING SEQUENCE:
!       Error_Status = CRTM_Atmosphere_ClearSkyCopy( Atm, Atm_Clear )
!
! INPUTS:
!       Atm:             Atmosphere structure to copy
!                        UNITS:      N/A
!                        TYPE:       CRTM_Atmosphere_type
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(IN)
!
! OUTPUTS:
!       Atm_Clear:       Copy of the input atmosphere but withut cloud information.
!                        UNITS:      N/A
!                        TYPE:       CRTM_Atmosphere_type
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(OUT)
!
! FUNCTION RESULT:
!       Error_Status:    The return value is an integer defining the error status.
!                        The error codes are defined in the Message_Handler module.
!                        If == SUCCESS the operation was successful
!                           == FAILURE an error occurred
!                        UNITS:      N/A
!                        TYPE:       INTEGER
!                        DIMENSION:  Scalar
!
!:sdoc-:
!--------------------------------------------------------------------------------

  FUNCTION CRTM_Atmosphere_ClearSkyCopy( atm, atm_clear ) RESULT( err_stat )
    ! Arguments
    TYPE(CRTM_Atmosphere_type), INTENT(IN)  :: atm
    TYPE(CRTM_Atmosphere_type), INTENT(OUT) :: atm_clear
    ! Function result
    INTEGER :: err_stat
    ! Local parameters
    CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'CRTM_Atmosphere_ClearSkyCopy'
    ! Local variables
    CHARACTER(ML) :: err_msg
    INTEGER :: i, k

    ! Set up
    err_stat = SUCCESS
    ! ...Check input
    IF ( .NOT. CRTM_Atmosphere_Associated(atm) ) THEN
      err_stat = FAILURE
      err_msg = 'Input Atmosphere structure not allocated'
      CALL Display_Message( ROUTINE_NAME, err_msg, err_stat )
      RETURN
    END IF

    ! Create the output structure
    CALL CRTM_Atmosphere_Create( atm_clear      , &
                                 atm%n_Layers   , &
                                 atm%n_Absorbers, &
                                 0              , &  ! NO CLOUDS !
                                 atm%n_Aerosols   )
    IF ( .NOT. CRTM_Atmosphere_Associated(atm_clear) ) THEN
      err_stat = FAILURE
      err_msg = 'Error allocating output Clear-Sky Atmosphere structure'
      CALL Display_Message( ROUTINE_NAME, err_msg, err_stat )
      RETURN
    END IF

    ! Copy over data
    ! ...Extra dimensions
    atm_clear%n_Added_Layers = atm%n_Added_Layers
    ! ...Layer independent data
    atm_clear%Climatology    = atm%Climatology
    atm_clear%Absorber_ID    = atm%Absorber_ID
    atm_clear%Absorber_Units = atm%Absorber_Units
    ! ...Layer dependent data
    k = atm%n_Layers
    atm_clear%Level_Pressure = atm%Level_Pressure(0:k)
    atm_clear%Pressure       = atm%Pressure(1:k)
    atm_clear%Temperature    = atm%Temperature(1:k)
    atm_clear%Absorber       = atm%Absorber(1:k,:)
    atm_clear%Cloud_Fraction = atm%Cloud_Fraction(1:k)
    ! ...Aerosol components
    IF ( atm%n_Aerosols > 0 ) THEN
      DO i = 1, atm%n_Aerosols
        atm_clear%Aerosol(i) = atm%Aerosol(i)
      END DO
    END IF

  END FUNCTION CRTM_Atmosphere_ClearSkyCopy

!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       CRTM_Atmosphere_ClearSkyCopy_TL
!
! PURPOSE:
!       Function to copy an instance of a tangent-linear CRTM Atmosphere object
!       but without the clouds included.
!!
! CALLING SEQUENCE:
!       Error_Status = CRTM_Atmosphere_ClearSkyCopy_TL( Atm, Atm_TL, Atm_Clear_TL )
!
! INPUTS:
!       Atm:             Atmosphere object for consistency checking
!                        UNITS:      N/A
!                        TYPE:       CRTM_Atmosphere_type
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(IN)
!
!       Atm_TL:          Tangent-linear Atmosphere object to copy. This object
!                        must be the tangent-linear equivalent of the input
!                        forward Atm object.
!                        This
!                        UNITS:      N/A
!                        TYPE:       CRTM_Atmosphere_type
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(IN)
!
! OUTPUTS:
!       Atm_Clear_TL:    Copy of the input atmosphere but withut cloud information.
!                        UNITS:      N/A
!                        TYPE:       CRTM_Atmosphere_type
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(OUT)
!
! FUNCTION RESULT:
!       Error_Status:    The return value is an integer defining the error status.
!                        The error codes are defined in the Message_Handler module.
!                        If == SUCCESS the operation was successful
!                           == FAILURE an error occurred
!                        UNITS:      N/A
!                        TYPE:       INTEGER
!                        DIMENSION:  Scalar
!
!:sdoc-:
!--------------------------------------------------------------------------------

  FUNCTION CRTM_Atmosphere_ClearSkyCopy_TL( &
    atm         , &  ! FWD input
    atm_TL      , &  ! TL  input
    atm_clear_TL) &  ! TL  output
  RESULT( err_stat )
    ! Arguments
    TYPE(CRTM_Atmosphere_type), INTENT(IN)  :: atm
    TYPE(CRTM_Atmosphere_type), INTENT(IN)  :: atm_TL
    TYPE(CRTM_Atmosphere_type), INTENT(OUT) :: atm_clear_TL
    ! Function result
    INTEGER :: err_stat
    ! Local parameters
    CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'CRTM_Atmosphere_ClearSkyCopy_TL'
    ! Local variables
    CHARACTER(ML) :: err_msg
    INTEGER :: i, k

    ! Set up
    err_stat = SUCCESS
    ! ...Check input allocation
    IF ( .NOT. CRTM_Atmosphere_Associated(atm   ) .OR. &
         .NOT. CRTM_Atmosphere_Associated(atm_TL) ) THEN
      err_stat = FAILURE
      err_msg = 'Input Atmosphere structures not allocated'
      CALL Display_Message( ROUTINE_NAME, err_msg, err_stat )
      RETURN
    END IF
    ! ...Dimension consistency
    IF ( (atm%n_Layers       /= atm_TL%n_Layers      ) .OR. &
         (atm%n_Absorbers    /= atm_TL%n_Absorbers   ) .OR. &
         (atm%n_Clouds       /= atm_TL%n_Clouds      ) .OR. &
         (atm%n_Aerosols     /= atm_TL%n_Aerosols    ) .OR. &
         (atm%n_Added_Layers /= atm_TL%n_Added_Layers) ) THEN
      err_stat = FAILURE
      err_msg = 'Input Atmosphere structures have incongruent dimensions'
      CALL Display_Message( ROUTINE_NAME, err_msg, err_stat )
      RETURN
    END IF
    ! ...Non-layer dependent data consistency
    IF (    (atm%Climatology    /= atm_TL%Climatology   ) .OR. &
         ANY(atm%Absorber_ID    /= atm_TL%Absorber_ID   ) .OR. &
         ANY(atm%Absorber_Units /= atm_TL%Absorber_Units) ) THEN
      err_stat = FAILURE
      err_msg = 'Input Atmosphere structures have incongruent layer independent data'
      CALL Display_Message( ROUTINE_NAME, err_msg, err_stat )
      RETURN
    END IF

    ! Create the output structure
    CALL CRTM_Atmosphere_Create( atm_clear_TL      , &
                                 atm_TL%n_Layers   , &
                                 atm_TL%n_Absorbers, &
                                 0                 , &  ! NO CLOUDS !
                                 atm_TL%n_Aerosols   )
    IF ( .NOT. CRTM_Atmosphere_Associated(atm_clear_TL) ) THEN
      err_stat = FAILURE
      err_msg = 'Error allocating output Clear-Sky Atmosphere structure'
      CALL Display_Message( ROUTINE_NAME, err_msg, err_stat )
      RETURN
    END IF

    ! Copy over data
    ! ...Extra dimensions
    atm_clear_TL%n_Added_Layers = atm_TL%n_Added_Layers
    ! ...Layer independent data
    atm_clear_TL%Climatology    = atm_TL%Climatology
    atm_clear_TL%Absorber_ID    = atm_TL%Absorber_ID
    atm_clear_TL%Absorber_Units = atm_TL%Absorber_Units
    ! ...Layer dependent data
    k = atm%n_Layers
    atm_clear_TL%Level_Pressure = atm_TL%Level_Pressure(0:k)
    atm_clear_TL%Pressure       = atm_TL%Pressure(1:k)
    atm_clear_TL%Temperature    = atm_TL%Temperature(1:k)
    atm_clear_TL%Absorber       = atm_TL%Absorber(1:k,:)
    atm_clear_TL%Cloud_Fraction = atm_TL%Cloud_Fraction(1:k)
    ! ...Aerosol components
    IF ( atm_TL%n_Aerosols > 0 ) THEN
      DO i = 1, atm_TL%n_Aerosols
        atm_clear_TL%Aerosol(i) = atm_TL%Aerosol(i)
      END DO
    END IF

  END FUNCTION CRTM_Atmosphere_ClearSkyCopy_TL

!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       CRTM_Atmosphere_ClearSkyCopy_AD
!
! PURPOSE:
!       Function to perform the adjoint copy of an instance of the CRTM
!       Atmosphere object without the clouds included.
!
! CALLING SEQUENCE:
!       Error_Status = CRTM_Atmosphere_ClearSkyCopy_AD( Atm, Atm_Clear_AD, Atm_AD )
!
! INPUTS:
!       Atm:             Atmosphere object for consistency checking
!                        UNITS:      N/A
!                        TYPE:       CRTM_Atmosphere_type
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(IN)
!
!       Atm_Clear_AD:    Adjoint Clear-Sky Atmosphere structure to copy
!                        UNITS:      N/A
!                        TYPE:       CRTM_Atmosphere_type
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(IN)
!
! OUTPUTS:
!       Atm_AD:          Adjoint copy of the input atmosphere. This object
!                        must be the adjoint equivalent of the input
!                        forward Atm object.
!                        UNITS:      N/A
!                        TYPE:       CRTM_Atmosphere_type
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(OUT)
!
! FUNCTION RESULT:
!       Error_Status:    The return value is an integer defining the error status.
!                        The error codes are defined in the Message_Handler module.
!                        If == SUCCESS the operation was successful
!                           == FAILURE an error occurred
!                        UNITS:      N/A
!                        TYPE:       INTEGER
!                        DIMENSION:  Scalar
!
!:sdoc-:
!--------------------------------------------------------------------------------

  FUNCTION CRTM_Atmosphere_ClearSkyCopy_AD( &
    atm         , &  ! FWD input
    atm_clear_AD, &  ! AD  input
    atm_AD      ) &  ! AD  output
  RESULT( err_stat )
    ! Arguments
    TYPE(CRTM_Atmosphere_type), INTENT(IN)     :: atm
    TYPE(CRTM_Atmosphere_type), INTENT(IN OUT) :: atm_clear_AD
    TYPE(CRTM_Atmosphere_type), INTENT(IN OUT) :: atm_AD
    ! Function result
    INTEGER :: err_stat
    ! Local parameters
    CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'CRTM_Atmosphere_ClearSkyCopy_AD'
    ! Local variables
    CHARACTER(ML) :: err_msg
    INTEGER :: i, k

    ! Set up
    err_stat = SUCCESS
    ! ...Check input allocation
    IF ( .NOT. CRTM_Atmosphere_Associated(atm         ) .OR. &
         .NOT. CRTM_Atmosphere_Associated(atm_clear_AD) .OR. &
         .NOT. CRTM_Atmosphere_Associated(atm_AD      ) ) THEN
      err_stat = FAILURE
      err_msg = 'Input Atmosphere structures not allocated'
      CALL Display_Message( ROUTINE_NAME, err_msg, err_stat )
      RETURN
    END IF
    ! ...Dimensional consistency
    IF ( (atm%n_Layers       /= atm_AD%n_Layers      ) .OR. &
         (atm%n_Absorbers    /= atm_AD%n_Absorbers   ) .OR. &
         (atm%n_Clouds       /= atm_AD%n_Clouds      ) .OR. &
         (atm%n_Aerosols     /= atm_AD%n_Aerosols    ) .OR. &
         (atm%n_Added_Layers /= atm_AD%n_Added_Layers) ) THEN
      err_stat = FAILURE
      err_msg = 'Input Atm and Atm_AD structures have incongruent dimensions'
      CALL Display_Message( ROUTINE_NAME, err_msg, err_stat )
      RETURN
    END IF
    IF ( (atm_clear_AD%n_Layers       /= atm_AD%n_Layers      ) .OR. &
         (atm_clear_AD%n_Absorbers    /= atm_AD%n_Absorbers   ) .OR. &
         (atm_clear_AD%n_Aerosols     /= atm_AD%n_Aerosols    ) .OR. &
         (atm_clear_AD%n_Clouds       /= 0                    ) .OR. &  ! NO CLOUDS !
         (atm_clear_AD%n_Added_Layers /= atm_AD%n_Added_Layers) ) THEN
      err_stat = FAILURE
      err_msg = 'Input Atm_Clear_AD structures has incongruent dimensions'
      CALL Display_Message( ROUTINE_NAME, err_msg, err_stat )
      RETURN
    END IF
    ! ...Non-layer dependent data consistency
    IF (    (atm%Climatology    /= atm_AD%Climatology   ) .OR. &
         ANY(atm%Absorber_ID    /= atm_AD%Absorber_ID   ) .OR. &
         ANY(atm%Absorber_Units /= atm_AD%Absorber_Units) .OR. &
            (atm%Climatology    /= atm_clear_AD%Climatology   ) .OR. &
         ANY(atm%Absorber_ID    /= atm_clear_AD%Absorber_ID   ) .OR. &
         ANY(atm%Absorber_Units /= atm_clear_AD%Absorber_Units) ) THEN
      err_stat = FAILURE
      err_msg = 'Atmosphere structures have incongruent layer independent data'
      CALL Display_Message( ROUTINE_NAME, err_msg, err_stat )
      RETURN
    END IF

    ! Adjoint copy of data
    ! ...Aerosol components
    IF ( atm%n_Aerosols > 0 ) THEN
      DO i = 1, atm%n_Aerosols
        atm_AD%Aerosol(i) = atm_AD%Aerosol(i) + atm_clear_AD%Aerosol(i)
      END DO
    END IF
    ! ...Layer dependent data
    k = atm%n_Layers
    atm_AD%Level_Pressure(0:k) = atm_AD%Level_Pressure(0:k) + atm_clear_AD%Level_Pressure(0:k)
    atm_AD%Pressure(1:k)       = atm_AD%Pressure(1:k)       + atm_clear_AD%Pressure(1:k)
    atm_AD%Temperature(1:k)    = atm_AD%Temperature(1:k)    + atm_clear_AD%Temperature(1:k)
    atm_AD%Absorber(1:k,:)     = atm_AD%Absorber(1:k,:)     + atm_clear_AD%Absorber(1:k,:)
    atm_AD%Cloud_Fraction(1:k) = atm_AD%Cloud_Fraction(1:k) + atm_clear_AD%Cloud_Fraction(1:k)

    ! Zero the clear result, as it has no more impact
    CALL CRTM_Atmosphere_Zero( atm_clear_AD )

  END FUNCTION CRTM_Atmosphere_ClearSkyCopy_AD

!##################################################################################
!##################################################################################
!##                                                                              ##
!##                          ## PRIVATE MODULE ROUTINES ##                       ##
!##                                                                              ##
!##################################################################################
!##################################################################################

  ! Subprogram to determine the number of extra layers required
  FUNCTION Extra_Layers( Atm ) RESULT( n )
    TYPE(CRTM_Atmosphere_type), INTENT(IN) :: Atm
    INTEGER :: n
    DO n = 1, SIZE(MODEL_LEVEL_PRESSURE)
      IF ( MODEL_LEVEL_PRESSURE(n) >= Atm%Level_Pressure(0) ) RETURN
    END DO
    n = 0
  END FUNCTION Extra_Layers

  ! Subprogram to compute the average layer pressure
  SUBROUTINE Layer_P( p, p_layer )
    REAL(fp), INTENT(IN)  :: p(2)     ! Input
    REAL(fp), INTENT(OUT) :: p_layer  ! Output
    p_layer = (p(2)-p(1))/LOG(p(2)/p(1))
  END SUBROUTINE Layer_P

  ! Subprogram to compute the average layer amount of X
  SUBROUTINE Layer_X( x, x_layer )
    REAL(fp), INTENT(IN)  :: x(2)
    REAL(fp), INTENT(OUT) :: x_layer
    x_layer = POINT_5*(x(1)+x(2))
  END SUBROUTINE Layer_X

  ! Subprogram to compute the interpolating polynomial linear in log(p)
  SUBROUTINE Interp_LPoly( p_int, p, lpoly )
    REAL(fp), INTENT(IN)  :: p_int
    REAL(fp), INTENT(IN)  :: p(2)
    REAL(fp), INTENT(OUT) :: lpoly
    lpoly = (LOG(p_int)-LOG(p(1))) / (LOG(p(2))-LOG(p(1)))
  END SUBROUTINE Interp_LPoly

  ! Subprogram to perform linear interpolation
  SUBROUTINE Interp_Linear( lpoly, x, x_int )
    REAL(fp), INTENT(IN)  :: lpoly
    REAL(fp), INTENT(IN)  :: x(2)
    REAL(fp), INTENT(OUT) :: x_int
    x_int = (x(2)-x(1))*lpoly + x(1)
  END SUBROUTINE Interp_Linear

  ! Subprogram to shifted the added profile layers
  SUBROUTINE Shift_Profile( lpoly, x_toa, x_shifted )
    REAL(fp), INTENT(IN)     :: lpoly
    REAL(fp), INTENT(IN)     :: x_toa(2)
    REAL(fp), INTENT(IN OUT) :: x_shifted(:)
    INTEGER :: n
    REAL(fp) :: x_int, dx
    n = SIZE(x_shifted)
    CALL Interp_Linear( lpoly, x_toa, x_int )
    dx = x_int - x_shifted(n)
    x_shifted = x_shifted + dx
  END SUBROUTINE Shift_Profile

END MODULE CRTM_Atmosphere

max_stars_repo_name: hsbadr/crtm

max_issues_repo_name: hsbadr/crtm

max_forks_repo_name: hsbadr/crtm


text:
!-------------------------------------------------------------------------
!Author: Matias Pablo, Borghi Orue
!Title: 2D Monte Carlo simulation of the Ising Model
!Version: 1
!Comments:  $ gfortran ranmar_seed.f95 subroutines.f95 ising.f95 -o runv1
!INPUTS: None
!OUTPUTS: Two dat files: 'Magnetization30L.dat' and 'Energy30L.dat'
!-------------------------------------------------------------------------
program ising

use subroutines
use rm
use iso_fortran_env
implicit none
	integer, parameter :: dp = selected_real_kind(15, 307)
	integer::ncols, average,iniseed
	integer::posx,posy,dim
	integer,allocatable::ip(:),im(:)
	integer:: MCS, EquilSteps
	integer:: i,j,k,indexTemp,printInfo !Dummy indeces
	integer::index1,index2
	real::Delta_E
	integer,allocatable::matriz(:,:)
	integer::promediateCounter = 0
	real::rvec(1)
	real(dp),allocatable::MagnetizationAver(:), MagnetizationAver2(:),MagnetizationAver4(:),EnergyAver(:),EnergyAver2(:)
	real,dimension(20)::DeltaE
	real,allocatable::alpha(:,:)
	real::alpha11,alpha33
	!Temp quantities
	real::beta
	!Field quantities
	real::currentField

	integer::fieldPoints,indexField
	character(len=10):: stringL
	character(len=100)::MagFile ='Magnetization.dat', EnergFile = 'Energy20L.dat'
	character(len=100):: folderName,FileOutput='OutputFile.log'

	real:: fieldValues1(1:40) = (/(i*0.002475, i=1,40,1) /),fieldValues2(1:40) = (/(i*0.002475+0.1,i=1,40,1) /)!,fieldValues3(-20:20) = (/(i*0.01,i=-20,20,1) /)
	real,allocatable::fieldValues(:)
	!real::fieldValues(1:30) = (/ (i*0.1,i=1,30,1) /)
	real:: currentTemp = 2.0,alpha13 = -1.0
	integer::nrows = 64

	!fieldPoints = size(fieldValues)

	! Read in input parameters from file "ising.in"
	open(unit=11,file="ConfigurationFile",status="old",action="read")
	read(11,*);read(11,*) printInfo
	read(11,*);read(11,*) iniseed
	read(11,*);read(11,*) MCS
	read(11,*);read(11,*) EquilSteps
	read(11,*);read(11,*) average
	read(11,*);read(11,*) alpha11
	read(11,*);read(11,*) alpha33
	close(11)

	!Fill the field values
	fieldPoints = size(fieldValues1) + size(fieldValues2) !+ size(fieldValues3)
	allocate(fieldValues(fieldPoints))
	fieldValues(1:size(fieldValues1)) = fieldValues1
	fieldValues(size(fieldValues1)+1:size(fieldValues1)+size(fieldValues2)) = fieldValues2
	!fieldValues(size(fieldValues1)+size(fieldValues2)+1:size(fieldValues)) = fieldValues3
	!print *, fieldValues
	!call EXIT(0)
	!-----------------
	ncols = nrows
	! Set the dimensions of the matrix of spin arrays. This program uses
	! periodic boundary conditions, so the first two rows and columns are
	! the same as the last two.
	allocate(matriz(nrows,ncols))
	allocate(alpha(3,3))

	allocate(MagnetizationAver(fieldPoints))
	allocate(MagnetizationAver2(fieldPoints))
	allocate(MagnetizationAver4(fieldPoints))
	allocate(EnergyAver(fieldPoints))
	allocate(EnergyAver2(fieldPoints))

	!Change this to take in count that the matrix can be rectangular
	allocate(im(nrows))
	allocate(ip(nrows))

	do i = 1,nrows
		ip(i) = i+1
		im(i) = i-1
	end do
	ip(nrows) = 1
	im(1) = nrows

	dim = nrows*ncols
	!MCS = 20000!dim*10
	!EquilSteps = 10000!int (0.5*MCS)

	write(stringL,'(I3.3)') nrows
	folderName = 'L='//trim(stringL)

	open(unit=2,file=trim(folderName)//'/'//trim(FileOutput),status="replace",action="write")
	write(2,*) 'Parameters in this simulation'
	write(2,*) 'L=',nrows
	write(2,*) 'Temp=',currentTemp
	write(2,*) 'Number of B points =', size(fieldValues)
	write(2,*) 'Alpha=',alpha11,' ',alpha13,' ',alpha33
	write(2,*) 'MCS=',MCS
	write(2,*) 'aver=',average

	if ( EquilSteps >= MCS ) then
		write(2,*) 'Error! Change the value of EquilSteps in the Configuration File! It can not be greater than MCS'
		call EXIT(0)
	end if

	write(2,*) 'Semilla RANMAR: ',iniseed
	write(2,*) "MCS = ",MCS, " EquilSteps = ", EquilSteps

	!If the file exists overwrite erase it.
	!open (UNIT=1, FILE=trim(folderName)//'/'//trim(MagFile),ACTION='write',STATUS='replace')
	!write(1,*) "Temp,B-Field, AverMagnetization, AverMagnetization^2, Suceptibility, Cumulant"
	!open (UNIT=33, FILE=trim(folderName)//'/'//trim(EnergFile),ACTION='write',STATUS='replace')
	!write(33,*) "Temp,B-Field, AverEnergy, AverEnergy^2, C_v"

	call ini_alpha(alpha,alpha11,alpha13,alpha33)
	!Initialize the matrix
	!with spin-up initial values
	call ini_matriz(matriz,nrows,ncols)
	close(2)

	if (printInfo /= 0) then
		call WriteSpinConf(matriz,nrows,ncols)
	end if

	!Initialize some variables to zero
	call SetToZeroVariables(MagnetizationAver,MagnetizationAver2,MagnetizationAver4,EnergyAver,EnergyAver2,fieldPoints)
	!call WriteSpinConf(matriz,nrows,ncols)
	beta = 1.0/currentTemp
	!Define the values of DeltaE knowing the Temperature
	!Only two possible values
	! 1) DeltaE(1) = exp(-4/Temp)
	! 2) DeltaE(2) = exp(-8/Temp)
	call Inicializacion_DeltaE(DeltaE,currentTemp)

	do k = 1 , average
		open(unit=2,file=trim(folderName)//'/'//trim(FileOutput),access="append",action="write")
		write(2,*) 'k = ',k,'/',average
		close(2)

		do indexField = 1,fieldPoints !Loop over all the temperatures from Tmax to Tmin
			promediateCounter = 0
			!Change the temperature value
			currentField = fieldValues(indexField)
			open(unit=2,file=trim(folderName)//'/'//trim(FileOutput),access="append",action="write")
			write(2,*), "Current Field B = ", currentField
			close(2)

			!Choose a random point (posx,posy) where we begin calculating
			call elegir_sitio(nrows,ncols,posx,posy,iniseed,k)

			!Sum over one MonteCarlo time //dim = 1MC
			!Sum over A MonteCarlo time: A*dim
			!In the next loop we flip a sping according to some conditions
			do i=1,MCS

				do j = 1,dim

					!Calculate DeltaE for a given point (posx,posy)
					call CalcDeltaEnergy(nrows,ncols,posx,posy,matriz,Delta_E,ip,im,alpha,currentField)

					!call ranmar(rvec,iniseed*verage)
					! ..Flipping Spin..
					!If the new Energy spin-state is negative with flip the spin
					!If not it will depend on the probability exp(-DeltaE/Temp) > random point
					if (Delta_E <= 0) then
						matriz(posx,posy) = - matriz(posx,posy)
					else
						call ranmar(rvec,1,iniseed*k)
						if (Delta_E == int(Delta_E) .AND. Delta_E<=20) then
							if (rvec(1) < DeltaE(int(Delta_E))) then
									matriz(posx,posy) = -matriz(posx,posy)
							end if
						else
							if ( rvec(1) < exp(-Delta_E/currentTemp) ) then
								matriz(posx,posy) = -matriz(posx,posy)
							end if
						end if

					end if

					!Choose a new random point (posx,posy) where we begin calculating
					call elegir_sitio(nrows,ncols,posx,posy,iniseed,k)

					!----------------------------------------------------------------
					!-----PRINTF IN DYNAMIC LATTICE ---------------------------------
					if (printInfo == 1) then
						call WriteSpinConf(matriz,nrows,ncols)
					end if
					!-----------------------------------------------------------------

				end do !End j loop

				!Calculate the Magnetization and Energy after equilibrium
				if (i > EquilSteps ) then
					promediateCounter = promediateCounter + 1
					call CalculationsINequilibrium (nrows,ncols,matriz,Delta_E,EnergyAver(indexField),EnergyAver2(indexField), &
										MagnetizationAver(indexField),MagnetizationAver2(indexField),MagnetizationAver4(indexField),ip,im,alpha,currentField)

					if (printInfo == 2) then
						call WriteSpinConf(matriz,nrows,ncols)
					else if (printInfo == 3 .AND. mod(i,MCS/10) == 0 .AND. k==1) then
						call WriteSpinConf(matriz,nrows,ncols)
					end if

				end if

			end do !End loop of the MC steps
			!call WriteSpinConf(matriz,nrows,ncols)
		end do ! End loop of Temperatures
	end do !End loop of averages
	!Write the results in a file
	call Write2File(promediateCounter, MagnetizationAver,MagnetizationAver2,MagnetizationAver4,EnergyAver,EnergyAver2,fieldValues,fieldPoints,average,currentTemp,folderName,MagFile)

	open(unit=2,file=trim(folderName)//'/'//trim(FileOutput),access="append",action="write")
	write(2,*), "Ended Program"
	close(2)

	!print *, 'promediateCounter = ',promediateCounter

	close(1)
	close(33)
end program ising

max_stars_repo_path: V6 Out-Of & On Equilibrium/V6.6 Magnetic Field/ising.f95

max_stars_repo_name: mattborghi/Thesis

max_issues_repo_name: mattborghi/Thesis

max_forks_repo_name: mattborghi/Thesis


text:
      SUBROUTINE PSPOSVX( FACT, UPLO, N, NRHS, A, IA, JA, DESCA, AF,
     $                    IAF, JAF, DESCAF, EQUED, SR, SC, B, IB, JB,
     $                    DESCB, X, IX, JX, DESCX, RCOND, FERR, BERR,
     $                    WORK, LWORK, IWORK, LIWORK, INFO )
*
*  -- ScaLAPACK routine (version 1.7) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     December 31, 1998
*
*     .. Scalar Arguments ..
      CHARACTER          EQUED, FACT, UPLO
      INTEGER            IA, IAF, IB, INFO, IX, JA, JAF, JB, JX, LIWORK,
     $                   LWORK, N, NRHS
      REAL               RCOND
*     ..
*     .. Array Arguments ..
      INTEGER            DESCA( * ), DESCAF( * ), DESCB( * ),
     $                   DESCX( * ), IWORK( * )
      REAL               A( * ), AF( * ),
     $                   B( * ), BERR( * ), FERR( * ),
     $                   SC( * ), SR( * ), WORK( * ), X( * )
*     ..
*
*  Purpose
*  =======
*
*  PSPOSVX uses the Cholesky factorization A = U**T*U or A = L*L**T to
*  compute the solution to a real system of linear equations
*
*        A(IA:IA+N-1,JA:JA+N-1) * X = B(IB:IB+N-1,JB:JB+NRHS-1),
*
*  where A(IA:IA+N-1,JA:JA+N-1) is an N-by-N matrix and X and
*  B(IB:IB+N-1,JB:JB+NRHS-1) are N-by-NRHS matrices.
*
*  Error bounds on the solution and a condition estimate are also
*  provided.  In the following comments Y denotes Y(IY:IY+M-1,JY:JY+K-1)
*  a M-by-K matrix where Y can be A, AF, B and X.
*
*  Notes
*  =====
*
*  Each global data object is described by an associated description
*  vector.  This vector stores the information required to establish
*  the mapping between an object element and its corresponding process
*  and memory location.
*
*  Let A be a generic term for any 2D block cyclicly distributed array.
*  Such a global array has an associated description vector DESCA.
*  In the following comments, the character _ should be read as
*  "of the global array".
*
*  NOTATION        STORED IN      EXPLANATION
*  --------------- -------------- --------------------------------------
*  DTYPE_A(global) DESCA( DTYPE_ )The descriptor type.  In this case,
*                                 DTYPE_A = 1.
*  CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
*                                 the BLACS process grid A is distribu-
*                                 ted over. The context itself is glo-
*                                 bal, but the handle (the integer
*                                 value) may vary.
*  M_A    (global) DESCA( M_ )    The number of rows in the global
*                                 array A.
*  N_A    (global) DESCA( N_ )    The number of columns in the global
*                                 array A.
*  MB_A   (global) DESCA( MB_ )   The blocking factor used to distribute
*                                 the rows of the array.
*  NB_A   (global) DESCA( NB_ )   The blocking factor used to distribute
*                                 the columns of the array.
*  RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
*                                 row of the array A is distributed.
*  CSRC_A (global) DESCA( CSRC_ ) The process column over which the
*                                 first column of the array A is
*                                 distributed.
*  LLD_A  (local)  DESCA( LLD_ )  The leading dimension of the local
*                                 array.  LLD_A >= MAX(1,LOCr(M_A)).
*
*  Let K be the number of rows or columns of a distributed matrix,
*  and assume that its process grid has dimension p x q.
*  LOCr( K ) denotes the number of elements of K that a process
*  would receive if K were distributed over the p processes of its
*  process column.
*  Similarly, LOCc( K ) denotes the number of elements of K that a
*  process would receive if K were distributed over the q processes of
*  its process row.
*  The values of LOCr() and LOCc() may be determined via a call to the
*  ScaLAPACK tool function, NUMROC:
*          LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
*          LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
*  An upper bound for these quantities may be computed by:
*          LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
*          LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
*
*  Description
*  ===========
*
*  The following steps are performed:
*
*  1. If FACT = 'E', real scaling factors are computed to equilibrate
*     the system:
*        diag(SR) * A * diag(SC) * inv(diag(SC)) * X = diag(SR) * B
*     Whether or not the system will be equilibrated depends on the
*     scaling of the matrix A, but if equilibration is used, A is
*     overwritten by diag(SR)*A*diag(SC) and B by diag(SR)*B.
*
*  2. If FACT = 'N' or 'E', the Cholesky decomposition is used to
*     factor the matrix A (after equilibration if FACT = 'E') as
*        A = U**T* U,  if UPLO = 'U', or
*        A = L * L**T,  if UPLO = 'L',
*     where U is an upper triangular matrix and L is a lower triangular
*     matrix.
*
*  3. The factored form of A is used to estimate the condition number
*     of the matrix A.  If the reciprocal of the condition number is
*     less than machine precision, steps 4-6 are skipped.
*
*  4. The system of equations is solved for X using the factored form
*     of A.
*
*  5. Iterative refinement is applied to improve the computed solution
*     matrix and calculate error bounds and backward error estimates
*     for it.
*
*  6. If equilibration was used, the matrix X is premultiplied by
*     diag(SR) so that it solves the original system before
*     equilibration.
*
*  Arguments
*  =========
*
*  FACT    (global input) CHARACTER
*          Specifies whether or not the factored form of the matrix A is
*          supplied on entry, and if not, whether the matrix A should be
*          equilibrated before it is factored.
*          = 'F':  On entry, AF contains the factored form of A.
*                  If EQUED = 'Y', the matrix A has been equilibrated
*                  with scaling factors given by S.  A and AF will not
*                  be modified.
*          = 'N':  The matrix A will be copied to AF and factored.
*          = 'E':  The matrix A will be equilibrated if necessary, then
*                  copied to AF and factored.
*
*  UPLO    (global input) CHARACTER
*          = 'U':  Upper triangle of A is stored;
*          = 'L':  Lower triangle of A is stored.
*
*  N       (global input) INTEGER
*          The number of rows and columns to be operated on, i.e. the
*          order of the distributed submatrix A(IA:IA+N-1,JA:JA+N-1).
*          N >= 0.
*
*  NRHS    (global input) INTEGER
*          The number of right hand sides, i.e., the number of columns
*          of the distributed submatrices B and X.  NRHS >= 0.
*
*  A       (local input/local output) REAL pointer into
*          the local memory to an array of local dimension
*          ( LLD_A, LOCc(JA+N-1) ).
*          On entry, the symmetric matrix A, except if FACT = 'F' and
*          EQUED = 'Y', then A must contain the equilibrated matrix
*          diag(SR)*A*diag(SC).  If UPLO = 'U', the leading
*          N-by-N upper triangular part of A contains the upper
*          triangular part of the matrix A, and the strictly lower
*          triangular part of A is not referenced.  If UPLO = 'L', the
*          leading N-by-N lower triangular part of A contains the lower
*          triangular part of the matrix A, and the strictly upper
*          triangular part of A is not referenced.  A is not modified if
*          FACT = 'F' or 'N', or if FACT = 'E' and EQUED = 'N' on exit.
*
*          On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by
*          diag(SR)*A*diag(SC).
*
*  IA      (global input) INTEGER
*          The row index in the global array A indicating the first
*          row of sub( A ).
*
*  JA      (global input) INTEGER
*          The column index in the global array A indicating the
*          first column of sub( A ).
*
*  DESCA   (global and local input) INTEGER array of dimension DLEN_.
*          The array descriptor for the distributed matrix A.
*
*  AF      (local input or local output) REAL pointer
*          into the local memory to an array of local dimension
*          ( LLD_AF, LOCc(JA+N-1)).
*          If FACT = 'F', then AF is an input argument and on entry
*          contains the triangular factor U or L from the Cholesky
*          factorization A = U**T*U or A = L*L**T, in the same storage
*          format as A.  If EQUED .ne. 'N', then AF is the factored form
*          of the equilibrated matrix diag(SR)*A*diag(SC).
*
*          If FACT = 'N', then AF is an output argument and on exit
*          returns the triangular factor U or L from the Cholesky
*          factorization A = U**T*U or A = L*L**T of the original
*          matrix A.
*
*          If FACT = 'E', then AF is an output argument and on exit
*          returns the triangular factor U or L from the Cholesky
*          factorization A = U**T*U or A = L*L**T of the equilibrated
*          matrix A (see the description of A for the form of the
*          equilibrated matrix).
*
*  IAF     (global input) INTEGER
*          The row index in the global array AF indicating the first
*          row of sub( AF ).
*
*  JAF     (global input) INTEGER
*          The column index in the global array AF indicating the
*          first column of sub( AF ).
*
*  DESCAF  (global and local input) INTEGER array of dimension DLEN_.
*          The array descriptor for the distributed matrix AF.
*
*  EQUED   (global input/global output) CHARACTER
*          Specifies the form of equilibration that was done.
*          = 'N':  No equilibration (always true if FACT = 'N').
*          = 'Y':  Equilibration was done, i.e., A has been replaced by
*                  diag(SR) * A * diag(SC).
*          EQUED is an input variable if FACT = 'F'; otherwise, it is an
*          output variable.
*
*  SR      (local input/local output) REAL array,
*                                                    dimension (LLD_A)
*          The scale factors for A distributed across process rows;
*          not accessed if EQUED = 'N'.  SR is an input variable if
*          FACT = 'F'; otherwise, SR is an output variable.
*          If FACT = 'F' and EQUED = 'Y', each element of SR must be
*          positive.
*
*  SC      (local input/local output) REAL array,
*                                              dimension (LOC(N_A))
*          The scale factors for A distributed across
*          process columns; not accessed if EQUED = 'N'. SC is an input
*          variable if FACT = 'F'; otherwise, SC is an output variable.
*          If FACT = 'F' and EQUED = 'Y', each element of SC must be
*          positive.
*
*  B       (local input/local output) REAL pointer into
*          the local memory to an array of local dimension
*          ( LLD_B, LOCc(JB+NRHS-1) ).
*          On entry, the N-by-NRHS right-hand side matrix B.
*          On exit, if EQUED = 'N', B is not modified; if TRANS = 'N'
*          and EQUED = 'R' or 'B', B is overwritten by diag(R)*B; if
*          TRANS = 'T' or 'C' and EQUED = 'C' or 'B', B is overwritten
*          by diag(C)*B.
*
*  IB      (global input) INTEGER
*          The row index in the global array B indicating the first
*          row of sub( B ).
*
*  JB      (global input) INTEGER
*          The column index in the global array B indicating the
*          first column of sub( B ).
*
*  DESCB   (global and local input) INTEGER array of dimension DLEN_.
*          The array descriptor for the distributed matrix B.
*
*  X       (local input/local output) REAL pointer into
*          the local memory to an array of local dimension
*          ( LLD_X, LOCc(JX+NRHS-1) ).
*          If INFO = 0, the N-by-NRHS solution matrix X to the original
*          system of equations.  Note that A and B are modified on exit
*          if EQUED .ne. 'N', and the solution to the equilibrated
*          system is inv(diag(SC))*X if TRANS = 'N' and EQUED = 'C' or
*          'B', or inv(diag(SR))*X if TRANS = 'T' or 'C' and EQUED = 'R'
*          or 'B'.
*
*  IX      (global input) INTEGER
*          The row index in the global array X indicating the first
*          row of sub( X ).
*
*  JX      (global input) INTEGER
*          The column index in the global array X indicating the
*          first column of sub( X ).
*
*  DESCX   (global and local input) INTEGER array of dimension DLEN_.
*          The array descriptor for the distributed matrix X.
*
*  RCOND   (global output) REAL
*          The estimate of the reciprocal condition number of the matrix
*          A after equilibration (if done).  If RCOND is less than the
*          machine precision (in particular, if RCOND = 0), the matrix
*          is singular to working precision.  This condition is
*          indicated by a return code of INFO > 0, and the solution and
*          error bounds are not computed.
*
*  FERR    (local output) REAL array, dimension (LOC(N_B))
*          The estimated forward error bounds for each solution vector
*          X(j) (the j-th column of the solution matrix X).
*          If XTRUE is the true solution, FERR(j) bounds the magnitude
*          of the largest entry in (X(j) - XTRUE) divided by
*          the magnitude of the largest entry in X(j).  The quality of
*          the error bound depends on the quality of the estimate of
*          norm(inv(A)) computed in the code; if the estimate of
*          norm(inv(A)) is accurate, the error bound is guaranteed.
*
*  BERR    (local output) REAL array, dimension (LOC(N_B))
*          The componentwise relative backward error of each solution
*          vector X(j) (i.e., the smallest relative change in
*          any entry of A or B that makes X(j) an exact solution).
*
*  WORK    (local workspace/local output) REAL array,
*                                                dimension (LWORK)
*          On exit, WORK(1) returns the minimal and optimal LWORK.
*
*  LWORK   (local or global input) INTEGER
*          The dimension of the array WORK.
*          LWORK is local input and must be at least
*          LWORK = MAX( PSPOCON( LWORK ), PSPORFS( LWORK ) )
*                  + LOCr( N_A ).
*          LWORK = 3*DESCA( LLD_ )
*
*          If LWORK = -1, then LWORK is global input and a workspace
*          query is assumed; the routine only calculates the minimum
*          and optimal size for all work arrays. Each of these
*          values is returned in the first entry of the corresponding
*          work array, and no error message is issued by PXERBLA.
*
*  IWORK   (local workspace/local output) INTEGER array,
*                                                  dimension (LIWORK)
*          On exit, IWORK(1) returns the minimal and optimal LIWORK.
*
*  LIWORK  (local or global input) INTEGER
*          The dimension of the array IWORK.
*          LIWORK is local input and must be at least
*          LIWORK = DESCA( LLD_ )
*          LIWORK = LOCr(N_A).
*
*          If LIWORK = -1, then LIWORK is global input and a workspace
*          query is assumed; the routine only calculates the minimum
*          and optimal size for all work arrays. Each of these
*          values is returned in the first entry of the corresponding
*          work array, and no error message is issued by PXERBLA.
*
*
*  INFO    (global output) INTEGER
*          = 0: successful exit
*          < 0: if INFO = -i, the i-th argument had an illegal value
*          > 0: if INFO = i, and i is
*               <= N: if INFO = i, the leading minor of order i of A
*                     is not positive definite, so the factorization
*                     could not be completed, and the solution and error
*                     bounds could not be computed.
*               = N+1: RCOND is less than machine precision.  The
*                     factorization has been completed, but the matrix
*                     is singular to working precision, and the solution
*                     and error bounds have not been computed.
*
*  =====================================================================
*
*     .. Parameters ..
      INTEGER            BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
     $                   LLD_, MB_, M_, NB_, N_, RSRC_
      PARAMETER          ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1,
     $                     CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6,
     $                     RSRC_ = 7, CSRC_ = 8, LLD_ = 9 )
      REAL               ONE, ZERO
      PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            EQUIL, LQUERY, NOFACT, RCEQU
      INTEGER            I, IACOL, IAROW, IAFROW, IBROW, IBCOL, ICOFF,
     $                   ICOFFA, ICTXT, IDUMM, IIA, IIB, IIX, INFEQU,
     $                   IROFF, IROFFA, IROFFAF, IROFFB, IROFFX, IXCOL,
     $                   IXROW, J, JJA, JJB, JJX, LDB, LDX, LIWMIN,
     $                   LWMIN, MYCOL, MYROW, NP, NPCOL, NPROW, NRHSQ,
     $                   NQ
      REAL               AMAX, ANORM, BIGNUM, SCOND, SMAX, SMIN, SMLNUM
*     ..
*     .. Local Arrays ..
      INTEGER            IDUM1( 5 ), IDUM2( 5 )
*     ..
*     .. External Subroutines ..
      EXTERNAL           BLACS_GRIDINFO, CHK1MAT, PCHK2MAT, INFOG2L,
     $                   PSPOCON, PSPOEQU, PSPORFS,
     $                   PSPOTRF, PSPOTRS,
     $                   PSLACPY, PSLAQSY, PXERBLA,
     $                   SGAMN2D, SGAMX2D
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      INTEGER            INDXG2P, NUMROC
      REAL               PSLANSY, PSLAMCH
      EXTERNAL           INDXG2P, LSAME, NUMROC, PSLANSY, PSLAMCH
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          ICHAR, MAX, MIN, MOD
*     ..
*     .. Executable Statements ..
*
*     Get grid parameters
*
      ICTXT = DESCA( CTXT_ )
      CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL )
*
*     Test the input parameters
*
      INFO = 0
      IF( NPROW.EQ.-1 ) THEN
         INFO = -(800+CTXT_)
      ELSE
         CALL CHK1MAT( N, 3, N, 3, IA, JA, DESCA, 8, INFO )
         IF( LSAME( FACT, 'F' ) )
     $      CALL CHK1MAT( N, 3, N, 3, IAF, JAF, DESCAF, 12, INFO )
         CALL CHK1MAT( N, 3, NRHS, 4, IB, JB, DESCB, 20, INFO )
         IF( INFO.EQ.0 ) THEN
            IAROW = INDXG2P( IA, DESCA( MB_ ), MYROW, DESCA( RSRC_ ),
     $                       NPROW )
            IAFROW = INDXG2P( IAF, DESCAF( MB_ ), MYROW,
     $                        DESCAF( RSRC_ ), NPROW )
            IBROW = INDXG2P( IB, DESCB( MB_ ), MYROW, DESCB( RSRC_ ),
     $                       NPROW )
            IXROW = INDXG2P( IX, DESCX( MB_ ), MYROW, DESCX( RSRC_ ),
     $                       NPROW )
            IROFFA = MOD( IA-1, DESCA( MB_ ) )
            IROFFAF = MOD( IAF-1, DESCAF( MB_ ) )
            ICOFFA = MOD( JA-1, DESCA( NB_ ) )
            IROFFB = MOD( IB-1, DESCB( MB_ ) )
            IROFFX = MOD( IX-1, DESCX( MB_ ) )
            CALL INFOG2L( IA, JA, DESCA, NPROW, NPCOL, MYROW,
     $                    MYCOL, IIA, JJA, IAROW, IACOL )
            NP = NUMROC( N+IROFFA, DESCA( MB_ ), MYROW, IAROW, NPROW )
            IF( MYROW.EQ.IAROW )
     $         NP = NP-IROFFA
            NQ = NUMROC( N+ICOFFA, DESCA( NB_ ), MYCOL, IACOL, NPCOL )
            IF( MYCOL.EQ.IACOL )
     $         NQ = NQ-ICOFFA
            LWMIN = 3*DESCA( LLD_ )
            LIWMIN = NP
            NOFACT = LSAME( FACT, 'N' )
            EQUIL = LSAME( FACT, 'E' )
            IF( NOFACT .OR. EQUIL ) THEN
               EQUED = 'N'
               RCEQU = .FALSE.
            ELSE
               RCEQU = LSAME( EQUED, 'Y' )
               SMLNUM = PSLAMCH( ICTXT, 'Safe minimum' )
               BIGNUM = ONE / SMLNUM
            END IF
            IF( .NOT.NOFACT .AND. .NOT.EQUIL .AND.
     $          .NOT.LSAME( FACT, 'F' ) ) THEN
               INFO = -1
            ELSE IF( .NOT.LSAME( UPLO, 'U' ) .AND.
     $               .NOT.LSAME( UPLO, 'L' ) ) THEN
               INFO = -2
            ELSE IF( IROFFA.NE.0 ) THEN
               INFO = -6
            ELSE IF( ICOFFA.NE.0 .OR. IROFFA.NE.ICOFFA ) THEN
               INFO = -7
            ELSE IF( DESCA( MB_ ).NE.DESCA( NB_ ) ) THEN
               INFO = -(800+NB_)
            ELSE IF( IAFROW.NE.IAROW ) THEN
               INFO = -10
            ELSE IF( IROFFAF.NE.0 ) THEN
               INFO = -10
            ELSE IF( ICTXT.NE.DESCAF( CTXT_ ) ) THEN
               INFO = -(1200+CTXT_)
            ELSE IF( LSAME( FACT, 'F' ) .AND. .NOT.
     $         ( RCEQU .OR. LSAME( EQUED, 'N' ) ) ) THEN
               INFO = -13
            ELSE
               IF( RCEQU ) THEN
*
                  SMIN = BIGNUM
                  SMAX = ZERO
                  DO 10 J = IIA, IIA + NP - 1
                     SMIN = MIN( SMIN, SR( J ) )
                     SMAX = MAX( SMAX, SR( J ) )
   10             CONTINUE
                  CALL SGAMN2D( ICTXT, 'Columnwise', ' ', 1, 1, SMIN,
     $                          1, IDUMM, IDUMM, -1, -1, MYCOL )
                  CALL SGAMX2D( ICTXT, 'Columnwise', ' ', 1, 1, SMAX,
     $                          1, IDUMM, IDUMM, -1, -1, MYCOL )
                  IF( SMIN.LE.ZERO ) THEN
                     INFO = -14
                  ELSE IF( N.GT.0 ) THEN
                     SCOND = MAX( SMIN, SMLNUM ) / MIN( SMAX, BIGNUM )
                  ELSE
                     SCOND = ONE
                  END IF
               END IF
            END IF
         END IF
*
         WORK( 1 ) = REAL( LWMIN )
         IWORK( 1 ) = LIWMIN
         LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
         IF( INFO.EQ.0 ) THEN
            IF( IBROW.NE.IAROW ) THEN
               INFO = -18
            ELSE IF( IXROW.NE.IBROW ) THEN
               INFO = -22
            ELSE IF( DESCB( MB_ ).NE.DESCA( NB_ ) ) THEN
               INFO = -(2000+NB_)
            ELSE IF( ICTXT.NE.DESCB( CTXT_ ) ) THEN
               INFO = -(2000+CTXT_)
            ELSE IF( ICTXT.NE.DESCX( CTXT_ ) ) THEN
               INFO = -(2400+CTXT_)
            ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
               INFO = -28
            ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
               INFO = -30
            END IF
            IDUM1( 1 ) = ICHAR( FACT )
            IDUM2( 1 ) = 1
            IDUM1( 2 ) = ICHAR( UPLO )
            IDUM2( 2 ) = 2
            IF( LSAME( FACT, 'F' ) ) THEN
               IDUM1( 3 ) = ICHAR( EQUED )
               IDUM2( 3 ) = 13
               IF( LWORK.EQ.-1 ) THEN
                  IDUM1( 4 ) = -1
               ELSE
                  IDUM1( 4 ) = 1
               END IF
               IDUM2( 4 ) = 28
               IF( LIWORK.EQ.-1 ) THEN
                  IDUM1( 5 ) = -1
               ELSE
                  IDUM1( 5 ) = 1
               END IF
               IDUM2( 5 ) = 30
               CALL PCHK2MAT( N, 3, N, 3, IA, JA, DESCA, 8, N, 3, NRHS,
     $                        4, IB, JB, DESCB, 19, 5, IDUM1, IDUM2,
     $                        INFO )
            ELSE
               IF( LWORK.EQ.-1 ) THEN
                  IDUM1( 3 ) = -1
               ELSE
                  IDUM1( 3 ) = 1
               END IF
               IDUM2( 3 ) = 28
               IF( LIWORK.EQ.-1 ) THEN
                  IDUM1( 4 ) = -1
               ELSE
                  IDUM1( 4 ) = 1
               END IF
               IDUM2( 4 ) = 30
               CALL PCHK2MAT( N, 3, N, 3, IA, JA, DESCA, 8, N, 3, NRHS,
     $                        4, IB, JB, DESCB, 19, 4, IDUM1, IDUM2,
     $                        INFO )
            END IF
         END IF
      END IF
*
      IF( INFO.NE.0 ) THEN
         CALL PXERBLA( ICTXT, 'PSPOSVX', -INFO )
         RETURN
      ELSE IF( LQUERY ) THEN
         RETURN
      END IF
*
      IF( EQUIL ) THEN
*
*        Compute row and column scalings to equilibrate the matrix A.
*
         CALL PSPOEQU( N, A, IA, JA, DESCA, SR, SC, SCOND, AMAX,
     $                 INFEQU )
*
         IF( INFEQU.EQ.0 ) THEN
*
*           Equilibrate the matrix.
*
            CALL PSLAQSY( UPLO, N, A, IA, JA, DESCA, SR, SC, SCOND,
     $                    AMAX, EQUED )
            RCEQU = LSAME( EQUED, 'Y' )
         END IF
      END IF
*
*     Scale the right-hand side.
*
      CALL INFOG2L( IB, JB, DESCB, NPROW, NPCOL, MYROW, MYCOL, IIB,
     $              JJB, IBROW, IBCOL )
      LDB =  DESCB( LLD_ )
      IROFF = MOD( IB-1, DESCB( MB_ ) )
      ICOFF = MOD( JB-1, DESCB( NB_ ) )
      NP = NUMROC( N+IROFF, DESCB( MB_ ), MYROW, IBROW, NPROW )
      NRHSQ = NUMROC( NRHS+ICOFF, DESCB( NB_ ), MYCOL, IBCOL, NPCOL )
      IF( MYROW.EQ.IBROW ) NP = NP-IROFF
      IF( MYCOL.EQ.IBCOL ) NRHSQ = NRHSQ-ICOFF
*
      IF( RCEQU ) THEN
         DO 30 J = JJB, JJB+NRHSQ-1
            DO 20 I = IIB, IIB+NP-1
               B( I + ( J-1 )*LDB ) = SR( I )*B( I + ( J-1 )*LDB )
   20       CONTINUE
   30    CONTINUE
      END IF
*
      IF( NOFACT .OR. EQUIL ) THEN
*
*        Compute the Cholesky factorization A = U'*U or A = L*L'.
*
         CALL PSLACPY( 'Full', N, N, A, IA, JA, DESCA, AF, IAF, JAF,
     $                 DESCAF )
         CALL PSPOTRF( UPLO, N, AF, IAF, JAF, DESCAF, INFO )
*
*        Return if INFO is non-zero.
*
         IF( INFO.NE.0 ) THEN
            IF( INFO.GT.0 )
     $         RCOND = ZERO
            RETURN
         END IF
      END IF
*
*     Compute the norm of the matrix A.
*
      ANORM = PSLANSY( '1', UPLO, N, A, IA, JA, DESCA, WORK )
*
*     Compute the reciprocal of the condition number of A.
*
      CALL PSPOCON( UPLO, N, AF, IAF, JAF, DESCAF, ANORM, RCOND, WORK,
     $              LWORK, IWORK, LIWORK, INFO )
*
*     Return if the matrix is singular to working precision.
*
      IF( RCOND.LT.PSLAMCH( ICTXT, 'Epsilon' ) ) THEN
         INFO = IA + N
         RETURN
      END IF
*
*     Compute the solution matrix X.
*
      CALL PSLACPY( 'Full', N, NRHS, B, IB, JB, DESCB, X, IX, JX,
     $              DESCX )
      CALL PSPOTRS( UPLO, N, NRHS, AF, IAF, JAF, DESCAF, X, IX, JX,
     $              DESCX, INFO )
*
*     Use iterative refinement to improve the computed solution and
*     compute error bounds and backward error estimates for it.
*
      CALL PSPORFS( UPLO, N, NRHS, A, IA, JA, DESCA, AF, IAF, JAF,
     $              DESCAF, B, IB, JB, DESCB, X, IX, JX, DESCX, FERR,
     $              BERR, WORK, LWORK, IWORK, LIWORK, INFO )
*
*     Transform the solution matrix X to a solution of the original
*     system.
*
      CALL INFOG2L( IX, JX, DESCX, NPROW, NPCOL, MYROW, MYCOL, IIX,
     $              JJX, IXROW, IXCOL )
      LDX = DESCX( LLD_ )
      IROFF = MOD( IX-1, DESCX( MB_ ) )
      ICOFF = MOD( JX-1, DESCX( NB_ ) )
      NP = NUMROC( N+IROFF, DESCX( MB_ ), MYROW, IXROW, NPROW )
      NRHSQ = NUMROC( NRHS+ICOFF, DESCX( NB_ ), MYCOL, IXCOL, NPCOL )
      IF( MYROW.EQ.IBROW ) NP = NP-IROFF
      IF( MYCOL.EQ.IBCOL ) NRHSQ = NRHSQ-ICOFF
*
      IF( RCEQU ) THEN
         DO 50 J = JJX, JJX+NRHSQ-1
            DO 40 I = IIX, IIX+NP-1
               X( I + ( J-1 )*LDX ) = SR( I )*X( I + ( J-1 )*LDX )
   40       CONTINUE
   50    CONTINUE
         DO 60 J = JJX, JJX+NRHSQ-1
            FERR( J ) = FERR( J ) / SCOND
   60    CONTINUE
      END IF
*
      WORK( 1 ) = REAL( LWMIN )
      IWORK( 1 ) = LIWMIN
      RETURN
*
*     End of PSPOSVX
*
      END

ext: f

max_stars_repo_name: jschueller/scalapack

max_issues_repo_name: jschueller/scalapack

max_forks_repo_name: jschueller/scalapack


text:
      SUBROUTINE SGENM (NTYPE,IFILE,SFILE,OFILE,ICODE,OCODE,CTYPES,
     1                  CTYPEO)
C
C     THIS SUBROUTINE MERGES CONVERTED SUBSTRUCTURING DATA WITH EXISTING
C     NASTRAN DATA
C
C     INPUTS
C     NTYPE  - NUMBER OF DIFFERENT SUBSTRUCTURING CARDS
C     IFILE  - INPUT FILE NAME
C     SFILE  - SCRATCH FILE NAME
C     OFILE  - OUTPUT FILE NAME
C     ICODE  - LOCATE CODES FOR INPUT CARD TYPES
C     OCODE  - LOCATE CODES FOR OUTPUT CARD TYPES
C     CTYPES - BCD NAMES OF SUBSTRUCTURING CARDS
C     CTYPEO - BCD NAMES OF CORRESPONDING NASTRAN CARDS
C
      INTEGER         SFILE,OFILE,ICODE(4,1),OCODE(4,1),CTYPES(2,8),
     1                CTYPEO(2,8),BUF1,BUF2,BUF3,Z,SYSBUF,OUTT,CARD(3),
     2                NLIMIT(3),SUBNAM(2)
      CHARACTER       UFM*23
      COMMON /XMSSG / UFM
      COMMON /BLANK / IDRY,NAME(2)
      COMMON /SGENCM/ NONO,NSS,IPTR,BUF1,BUF2,BUF3,NZ
      COMMON /ZZZZZZ/ Z(1)
      COMMON /SYSTEM/ SYSBUF,OUTT
      DATA    NLIMIT/ 3*2147483647 /
      DATA    SUBNAM/ 4HSGEN,4HM    /
C
C     OPEN FILES
C
      CALL GOPEN (IFILE,Z(BUF1),0)
      CALL GOPEN (SFILE,Z(BUF2),0)
      CALL GOPEN (OFILE,Z(BUF3),1)
C
C     READ HEADER FROM IFILE - DETERMINE IF SUBSTRUCTURING OR NASTRAN
C     CARD
C
      FILE = IFILE
   10 CALL READ (*1002,*1003,IFILE,CARD,3,0,IDX)
      IF (CARD(1) .EQ. NLIMIT(1)) GO TO 70
      DO 20 I = 1,NTYPE
      IF (ICODE(1,I).NE.CARD(1) .OR. ICODE(2,I).NE.CARD(2)) GO TO 20
C
C     SKIP RECORD IF SUBSTRUCTURING CARD
C
      CALL FWDREC (*70,IFILE)
      GO TO 10
   20 CONTINUE
      DO 30 I = 1,NTYPE
      IF (OCODE(1,I).NE.CARD(1) .OR. OCODE(2,I).NE.CARD(2)) GO TO 30
C
C     FATAL ERROR IF BOTH SUBSTRUCTURING AND NASTRAN CARDS
C
      IF (ICODE(4,I) .EQ. 0) GO TO 40
      NONO = 1
      J = OCODE(4,I)
      WRITE (OUTT,6330) UFM,NAME,(CTYPES(K,J),K=1,2),(CTYPEO(K,J),K=1,2)
      CALL FWDREC (*70,IFILE)
      GO TO 10
   30 CONTINUE
C
C     COPY RECORD FROM IFILE TO OUTPUT
C
   40 CALL WRITE (OFILE,CARD,3,0)
   50 CALL READ  (*1002,*60,IFILE,Z,NZ,0,NWDS)
      CALL WRITE (OFILE,Z,NZ,0)
      GO TO 50
   60 CALL WRITE (OFILE,Z,NWDS,1)
      GO TO 10
C
C     COPY RECORD FROM SFILE TO OUTPUT
C
   70 I1 = 1
   80 DO 90 I = I1,NTYPE
      IF (ICODE(4,I) .EQ. 1) GO TO 100
   90 CONTINUE
      GO TO 150
  100 CALL FREAD (SFILE,CARD,3,0)
      CALL WRITE (OFILE,OCODE(1,I),3,0)
      FILE = SFILE
  110 CALL READ  (*1002,*120,SFILE,Z,NZ,0,NWDS)
      CALL WRITE (OFILE,Z,NZ,0)
      GO TO 110
  120 CALL WRITE (OFILE,Z,NWDS,1)
      I1 = I + 1
      GO TO 80
C
C     CLOSE FILES
C
  150 CALL WRITE (OFILE,NLIMIT,3,1)
      CALL CLOSE (IFILE,1)
      CALL CLOSE (SFILE,1)
      CALL CLOSE (OFILE,1)
      RETURN
C
C     ERRORS
C
 1002 M = -2
      GO TO 2000
 1003 M = -3
 2000 CALL MESAGE (M,FILE,SUBNAM)
      RETURN
C
 6330 FORMAT (A23,' 6330, SOLUTION SUBSTRUCTURE ',2A4,3H - ,2A4,' AND ',
     1       2A4,' CARDS CANNOT BE USED TOGETHER.', /30X,
     2       'USE EITHER ONE, BUT NOT BOTH.')
      END


text:
program problem47
  use euler
  implicit none
  integer, parameter             :: window=4,maxprime=1000000
  integer*8                      :: i,j,k,distinct,nprimes
  integer*8, dimension(maxprime) :: primes

  nprimes=1
  do i=2,maxprime
     if (is_prime(i)) then
        primes(nprimes)=i
        nprimes=nprimes+1
     end if
  end do

  do i=1,huge(i) ! start of window of numbers
     do j=i,i+window-1 ! for each number in the window
        distinct=0
        do k=1,nprimes
           if (primes(k) > j) exit
           if (mod(j,primes(k)) == 0) distinct=distinct+1
           if (distinct > window) go to 1
        end do
        if (distinct < window) go to 1
     end do
     exit ! made it through all of the checks, answer found
1 end do

  print *, i

end program problem47

max_stars_repo_name: jamesmcclain/Dioskouroi

max_issues_repo_name: jamesmcclain/Dioskouroi

max_forks_repo_name: jamesmcclain/Dioskouroi


text:
In addition to being one of the owners of Watermelon Music, Jeff is a great dad and a very tasteful bass player.

20080127 10:13:23 nbsp Welcome to the Wiki. Your comment on Watermelon Music was really well written and if I had any musical talent at all Id shop there. The wiki community, of which you are a part, continues to struggle with how to handle reviews, both the negative and the suspiciously good. Your input along the way can help to make the wiki a better place as time goes by. Users/JasonAller

20090224 21:15:00 nbsp you tell em Jeff! Users/StevenDaubert

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


module procedures
implicit none
contains

subroutine writes2d (A)
   real(kind=8), allocatable, dimension(:,:), intent(in) :: A
   integer :: n, m, i, j

   n = size(A, 1)
   m = size(A, 2)

   do i = 1,n
      write(*,'(20F8.4)') (A(i,j), j = 1, m)
   end do
   write(*,*)
end subroutine writes2d

subroutine writes1d (A)
   real(kind=8), allocatable, dimension(:), intent(in) :: A
   integer :: n, i

   n = size(A)

   write(*,'(20F8.4)') (A(i), i = 1, n)
   write(*,*)
end subroutine writes1d

real(kind = 8) function MaxEig(A)
   real(kind=8), allocatable, dimension(:,:), intent(in) :: A
   real(kind=8), allocatable, dimension(:,:) :: Y, Y_old, Y_new
   real(kind=8) :: normalization_constant
   integer :: n, i, j

   !Initialize column vector:
   n = size(A,1)
   allocate (Y(n,1), Y_old(n,1))
   do i = 1, n
      Y(i,1) = 1.0d0
      Y_old(i,1) = 2.0d0 !initialize loop
   end do
   write(*,*) "initial vector:"
   call writes2d(Y)

   !perform calculation:
   i = 1
   do while (maxval(dabs(Y_old-Y)) .ge. 0.0001)
      do j = 1, n ! save old vector
         Y_old(j,1) = Y(j,1)
      end do
      Y = matmul(A,Y) !determine new vector
      MaxEig = maxval(Y)
      do j = 1, n
         Y(j,1) = Y(j,1)/MaxEig
      end do
      write(*,*) "iteration:", i
      write(*,'(A14,F10.6)') " tolerance:   ", maxval(dabs(Y_old-Y))
      write(*,'(A15,F8.4)') " eigenvalue:   ", MaxEig
      write(*,*) "eigenvector:"
      call writes2d(Y)
      i = i+1
   end do

   ! final normalization:
   normalization_constant = sqrt(Y(1,1)*Y(1,1) + Y(2,1)*Y(2,1) + Y(3,1)*Y(3,1))
   do i = 1, n
      Y(i,1) = Y(i,1)/normalization_constant
   end do
   write(*,*)
   write(*,'(A30)') " -----------------------------"
   write(*,'(A30)') "         final result         "
   write(*,'(A30)') " -----------------------------"
   write(*,*) "eigenvalue:", MaxEig
   write(*,*) "eigenvector, normalized to length 1 for comparison with matlab result:"
   call writes2d(Y)

   deallocate (Y, Y_old)
end function MaxEig

end module procedures


text: Their website claims that Toftware is a consulting firm providing custom software and database solutions through a wide variety of services.  We can take your project from concept to product, add new features to your current applications, or join your team to help bring an ongoing project to completion.  Bring us your vision, ideas, and goals.  We have the skills and experience to provide you with a low risk and high quality solution.


text:
      subroutine isource(sfcden,threshold,sx,sy,co2ice,fs)

!     Legacy Mars GCM v24
!     Mars Climate Modeling Center
!     NASA Ames Research Center

C  Dust tracer
C  Compute surface dust flux, in a fully interactive manner
C  interactive source
C  April 2002
C
C  SX         - surface stress x-direction  (N/m^2)
C  SY         - surface stress y-direction  (N/m^2)
C  CO2ICE     - Amount of CO2 ice on the ground
C  FS         - mass flux of dust lifted from the surface (kg m^-2 s^-1)
C  THRESHOLD  - stress threshold
C----------------------------------------------------------------------C

      use constants_h, only: GRAV
      implicit none

      real*8 sx, sy, co2ice, fs
      real*8 stress, sfcden ,threshold, hflux, alfa

!======================================================================!

c      alfa = 1.0e-5  ! Global dust storm
c      alfa = 3.0e-6  ! 24*16 dust interactive, no clouds
      alfa = 5.e-7  ! 24*16 dust interactive, interactive clouds

      stress = sqrt(sx**2 + sy**2)

C  No lifting if stress is below the threshold value (THRESHOLD), or if
C  there is CO2 ice on the ground (TG < 0).

      if(stress.le.threshold .or. co2ice.gt.0.0) then
        fs = 0.0
      else
C First the horizontal saltation flux
        hflux = 2.61 * ( stress**1.5/grav/sqrt(sfcden) ) *
     &           ( 1. - sqrt(threshold/stress) ) *
     &           ( 1. + sqrt(threshold/stress) ) ** 2
C Then the saltation flux
        fs = alfa*hflux
      end if

      return
      end

max_stars_repo_name: nasa/legacy-mars-global-climate-model

max_issues_repo_name: nasa/legacy-mars-global-climate-model

max_forks_repo_name: nasa/legacy-mars-global-climate-model


text:
!
! Copyright (C) 2001-2010 Quantum ESPRESSO group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!-----------------------------------------------------------------------
SUBROUTINE rgen ( dtau, rmax, mxr, at, bg, r, r2, nrm)
  !-----------------------------------------------------------------------
  !
  !   generates neighbours shells (cartesian, in units of lattice parameter)
  !   with length < rmax,and returns them in order of increasing length:
  !      r(:) = i*a1(:) + j*a2(:) + k*a3(:) - dtau(:),   r2 = r^2
  !   where a1, a2, a3 are primitive lattice vectors. Other input variables:
  !     mxr = maximum number of vectors
  !     at  = lattice vectors ( a1=at(:,1), a2=at(:,2), a3=at(:,3) )
  !     bg  = reciprocal lattice vectors ( b1=bg(:,1), b2=bg(:,2), b3=bg(:,3) )
  !   Other output variables:
  !     nrm = the number of vectors with r^2 < rmax^2
  !
  USE kinds, ONLY : DP
  !
  IMPLICIT NONE
  INTEGER, INTENT(in) :: mxr
  INTEGER, INTENT(out):: nrm
  REAL(DP), INTENT(in) :: at(3,3), bg(3,3), dtau(3), rmax
  REAL(DP), INTENT(out):: r(3,mxr), r2(mxr)
  !
  !    and here the local variables
  !
  INTEGER, ALLOCATABLE :: irr (:)
  INTEGER ::  nm1, nm2, nm3, i, j, k, ipol, ir, indsw, iswap
  real(DP) :: ds(3), dtau0(3)
  real(DP) :: t (3), tt, swap
  real(DP), EXTERNAL :: dnrm2
  !
  !
  nrm = 0
  IF (rmax==0.d0) RETURN

  ! bring dtau into the unit cell centered on the origin - prevents trouble
  ! if atomic positions are not centered around the origin but displaced
  ! far away (remember that translational invariance allows this!)
  !
  ds(:) = matmul( dtau(:), bg(:,:) )
  ds(:) = ds(:) - anint(ds(:))
  dtau0(:) = matmul( at(:,:), ds(:) )
  !
  ALLOCATE (irr( mxr))
  !
  ! these are estimates of the maximum values of needed integer indices
  !
  nm1 = int (dnrm2 (3, bg (1, 1), 1) * rmax) + 2
  nm2 = int (dnrm2 (3, bg (1, 2), 1) * rmax) + 2
  nm3 = int (dnrm2 (3, bg (1, 3), 1) * rmax) + 2
  !
  DO i = -nm1, nm1
     DO j = -nm2, nm2
        DO k = -nm3, nm3
           tt = 0.d0
           DO ipol = 1, 3
              t (ipol) = i*at (ipol, 1) + j*at (ipol, 2) + k*at (ipol, 3) &
                       - dtau0(ipol)
              tt = tt + t (ipol) * t (ipol)
           ENDDO
           IF (tt<=rmax**2.and.abs (tt) >1.d-10) THEN
              nrm = nrm + 1
              IF (nrm>mxr) CALL errore ('rgen', 'too many r-vectors', nrm)
              DO ipol = 1, 3
                 r (ipol, nrm) = t (ipol)
              ENDDO
              r2 (nrm) = tt
           ENDIF
        ENDDO
     ENDDO
  ENDDO
  !
  !   reorder the vectors in order of increasing magnitude
  !
  !   initialize the index inside sorting routine
  !
  irr (1) = 0
  IF (nrm>1) CALL hpsort (nrm, r2, irr)
  DO ir = 1, nrm - 1
20   indsw = irr (ir)
     IF (indsw/=ir) THEN
        DO ipol = 1, 3
           swap = r (ipol, indsw)
           r (ipol, indsw) = r (ipol, irr (indsw) )
           r (ipol, irr (indsw) ) = swap
        ENDDO
        iswap = irr (ir)
        irr (ir) = irr (indsw)
        irr (indsw) = iswap
        GOTO 20
     ENDIF

  ENDDO
  DEALLOCATE(irr)
  !
  RETURN
END SUBROUTINE rgen

max_stars_repo_name: utdsimmons/ohpc

max_issues_repo_name: utdsimmons/ohpc

max_forks_repo_name: utdsimmons/ohpc


text:
C File: test_module_11.f
C Purpose: Illustrates the use of multiple Fortran modules to define values
C    for several variables that are used in the program.  This program
C    differs from test_module_03.f in having modules that have "USE MODULE"s
C    inside them , as well as some variables that are declared to be PRIVATE.
C
C Compile and run this program as follows:
C
C    gfortran -c test_module_11.f   # << this creates the files "mymod11a.mod"
C                                                               "mymod11b.mod"
C                                                               "mymod11c.mod"
C    gfortran test_module_11.f      # << this creates a file "a.out"
C
C The output generated by this program is (see test_module_11-OUT.txt):
C
C mymod11B: X = 123; Y = 234; Z = 345; U = 567; V = 678
C mymod11C: X = 123; Y = 234; Z =   12.345; U = 567; V = 678
C pgm main: X = 123; Y =  98.765; Z =  12.345; U =  87.654; V = 678

!-------------------------------------------------------------------------------
      MODULE MYMOD11A
      IMPLICIT NONE

      INTEGER :: X = 123, Y = 234, Z = 345
      END MODULE mymod11a

!-------------------------------------------------------------------------------
      MODULE MYMOD11B
      USE MYMOD11A    ! << Note that "use mymod11a" is inside mymod11b

          private :: Z    ! << hides Z imported from mymod11a
          INTEGER :: U = 567, V = 678
      contains
          subroutine print_mymod11b()
 10       format ('mymod11B: X = ', I3, '; Y = ', I3, '; Z = ', I3,
     &                    '; U = ', I3, '; V = ', I3)
          write (*,10) X, Y, Z, U, V
          end subroutine print_mymod11b

      END MODULE mymod11b

!-------------------------------------------------------------------------------
      MODULE MYMOD11C
      USE mymod11B    ! << Note that "use mymod11b" is inside mymod11c

          PRIVATE :: Y, U ! << hides Y from mymod11a, U from mymod11b
          REAL :: Z = 12.345  ! << reusing Z
      contains
          subroutine print_mymod11c()
          call print_mymod11b()
 11       FORMAT('mymod11C: X = ', I3, '; Y = ', I3,'; Z = ', F8.3,
     &                     '; U = ', I3, '; V = ',I3)
          WRITE (*,11) X, Y, Z, U, V
          end subroutine print_mymod11c

      END MODULE mymod11C

!-------------------------------------------------------------------------------
      PROGRAM PGM
      USE mymod11c

      REAL :: Y = 98.765, U = 87.654  ! << reusing Y, U

      call print_mymod11c()

 10   FORMAT('pgm main: X = ', I3, '; Y = ', F7.3,'; Z = ', F7.3,
     &                     '; U = ', F7.3, '; V = ',I3)
      write (*,10) X, Y, Z, U, V

      stop
      end program PGM
!-------------------------------------------------------------------------------