#INCLUDE 'MR_H_ALIGN_PADDING.H'
!***********************************************************************************************************************************
! UNIT:
!
!  (MODULE)
!
! PURPOSE:
!
!   TO
!
! DEFINITION OF VARIABLES:
!
!
!
! RECORD OF REVISIONS:
!
!      DATE       |    PROGRAMMER    |    DESCRIPTION OF CHANGE
!      ====       |    ==========    |    =====================
!   20XX-XX-XX    |     DR. HYDE     |    ORIGINAL CODE.
!
!***********************************************************************************************************************************
  MODULE MR_MOD_GEN_INI_ZB

    USE MR_ERRORS_NONLINEAR_EQN_SOLVE

    USE MR_KINDS

    USE MR_MAC_PI

    USE MR_DEF_RANKS
    USE MR_DEF_CONSTS_N_REF_PARS
    USE MR_DEF_CURVED_GEOS

    IMPLICIT NONE

    PRIVATE

    PUBLIC :: MR_GEN_INI_ZB

    REAL   (GJRD_KIND) , ALLOCATABLE , DIMENSION(:) :: DAREA
    REAL   (GJRD_KIND) , ALLOCATABLE , DIMENSION(:) :: WTH , RTH
    REAL   (GJRD_KIND) , ALLOCATABLE , DIMENSION(:) :: RTH_BK_L , RTH_BK_R

    REAL   (PARD_KIND) , ALLOCATABLE , DIMENSION(:) :: DZB_BK_L , DZB_BK_R

    REAL   (CARD_KIND) , ALLOCATABLE , DIMENSION(:) :: ALFA_LBOUND , ALFA_RBOUND

    REAL   (CARD_KIND) , ALLOCATABLE , DIMENSION(:) :: ALFA , BETA , GAMA

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

  CONTAINS

!***********************************************************************************************************************************
! UNIT:
!
!  (SUBROUTINE)
!
! PURPOSE:
!
!   TO
!
! DEFINITION OF VARIABLES:
!
!
!
! RECORD OF REVISIONS:
!
!      DATE       |    PROGRAMMER    |    DESCRIPTION OF CHANGE
!      ====       |    ==========    |    =====================
!   20XX-XX-XX    |     DR. HYDE     |    ORIGINAL CODE.
!
!***********************************************************************************************************************************
  SUBROUTINE MR_INIT_WTH_N_RTH

    IMPLICIT NONE

    INTEGER(IJID_KIND) :: I , J

    ALLOCATE( DAREA(1:NI1(NI,GJRD_KIND)) )

      DAREA = 0
      DO J = 1 , NJ
        !DIR$ VECTOR ALIGNED
        DO I = 1 , NI
          DAREA( I ) = DAREA( I ) + MW( I , J )
        END DO
      END DO

      !DIR$ VECTOR ALIGNED
      DO I = 1 , NI
        WTH( I ) = 2.0 * DAREA( I )   &
        / ( SQRT( GVV( I , 0 ,1,1) ) + SQRT( GVV( I ,NJ ,1,1) ) )
      END DO

    DEALLOCATE( DAREA )

   !DIR$ VECTOR ALIGNED
    DO I = 1 , NI

      RTH( I ) =   &
      ( ( SIGN( ACOS( JOO(I-1, 0 ,1,1) / SQRT(GOO(I-1, 0 ,1,1)) ) , JOO(I-1, 0 ,2,1) )   &
        + SIGN( ACOS( JOO(I-1,NJ ,1,1) / SQRT(GOO(I-1,NJ ,1,1)) ) , JOO(I-1,NJ ,2,1) )   &
        )   &
      - ( SIGN( ACOS( JOO( I , 0 ,1,1) / SQRT(GOO( I , 0 ,1,1)) ) , JOO( I , 0 ,2,1) )   &
        + SIGN( ACOS( JOO( I ,NJ ,1,1) / SQRT(GOO( I ,NJ ,1,1)) ) , JOO( I ,NJ ,2,1) )   &
        )   &
      ) / ( SQRT( GVV( I , 0 ,1,1) ) + SQRT( GVV( I ,NJ ,1,1) ) )

      RTH( I ) = 1.0 / SIGN( MAX( ABS(RTH( I )) , 10.0*EPSILON(RTH( I )) ) , RTH( I ) )

      RTH_BK_L( I ) = RTH( I ) + 0.5 * WTH( I )
      RTH_BK_R( I ) = RTH( I ) - 0.5 * WTH( I )

    END DO

  END SUBROUTINE MR_INIT_WTH_N_RTH

!***********************************************************************************************************************************
! UNIT:
!
!  (SUBROUTINE)
!
! PURPOSE:
!
!   TO
!
! DEFINITION OF VARIABLES:
!
!
!
! RECORD OF REVISIONS:
!
!      DATE       |    PROGRAMMER    |    DESCRIPTION OF CHANGE
!      ====       |    ==========    |    =====================
!   20XX-XX-XX    |     DR. HYDE     |    ORIGINAL CODE.
!
!***********************************************************************************************************************************
  SUBROUTINE MR_GEN_INI_ZB( HTH , DZB_BK_MIN , DZB_BK_MAX , XI0 , XXIM , NBENDS , NI , NJ , ZB , ERROR , ERRMSG )

    IMPLICIT NONE

    REAL   (PARD_KIND) , INTENT(IN ) :: HTH
    REAL   (PARD_KIND) , INTENT(IN ) :: DZB_BK_MIN , DZB_BK_MAX
    REAL   (PARD_KIND) , INTENT(IN ) :: XI0 , XXIM

    INTEGER            , INTENT(IN ) :: NBENDS

    INTEGER(IJID_KIND) , INTENT(IN ) :: NI , NJ

    REAL   (FDRD_KIND) , INTENT(OUT) , DIMENSION(1:NI1(NI,FDRD_KIND),1:NJ) :: ZB

    INTEGER            , INTENT(OUT) :: ERROR
    CHARACTER(   *   ) , INTENT(OUT) :: ERRMSG

    REAL   (PARD_KIND) :: COEFFI_XXIM
    REAL   (PARD_KIND) :: COEFFI_DZB_BK_MIN_N_MAX_1 , COEFFI_DZB_BK_MIN_N_MAX_2

    REAL   (PARD_KIND) :: DYJ

    INTEGER(IJID_KIND) :: I , J

    COEFFI_XXIM = ( 0.25 - XXIM ) / ( 1.0 + SIN( 2.0 * PI * XXIM ) )

    COEFFI_DZB_BK_MIN_N_MAX_1 = 0.25 * ( DZB_BK_MAX + DZB_BK_MIN )
    COEFFI_DZB_BK_MIN_N_MAX_2 = 0.50 * ( DZB_BK_MAX - DZB_BK_MIN )

    IF( (.NOT. ALLOCATED(DZB_BK_L) ) .AND. (.NOT. ALLOCATED(DZB_BK_R) ) ) THEN
      ALLOCATE( DZB_BK_L(1:NI1(NI,PARD_KIND)) , DZB_BK_R(1:NI1(NI,PARD_KIND)) )
      CALL MR_GEN_INI_DZB_BK( XI0 , COEFFI_XXIM , COEFFI_DZB_BK_MIN_N_MAX_1 , COEFFI_DZB_BK_MIN_N_MAX_2 ,   &
      & NBENDS , NI , DZB_BK_L , DZB_BK_R , ERROR , ERRMSG )
      IF( ERROR < 0 ) THEN
        DEALLOCATE( DZB_BK_L , DZB_BK_R )
        RETURN
      END IF
    END IF

      IF( (.NOT. ALLOCATED(WTH) ) .AND. (.NOT. ALLOCATED(RTH) )   &
      &  .AND. (.NOT. ALLOCATED(RTH_BK_L) ) .AND. (.NOT. ALLOCATED(RTH_BK_R) ) ) THEN
        ALLOCATE( WTH(1:NI1(NI,GJRD_KIND)) , RTH(1:NI1(NI,GJRD_KIND)) ,   &
        & RTH_BK_L(1:NI1(NI,GJRD_KIND)) , RTH_BK_R(1:NI1(NI,GJRD_KIND)) )
        CALL MR_INIT_WTH_N_RTH
      END IF

        ALLOCATE( ALFA(1:NI1(NI,CARD_KIND)) , BETA(1:NI1(NI,CARD_KIND)) , GAMA(1:NI1(NI,CARD_KIND)) )

          IF( (.NOT. ALLOCATED(ALFA_LBOUND) ) .AND. (.NOT. ALLOCATED(ALFA_RBOUND) ) ) THEN
            ALLOCATE( ALFA_LBOUND(1:NI1(NI,CARD_KIND)) , ALFA_RBOUND(1:NI1(NI,CARD_KIND)) )
            ALFA_LBOUND = -8.001 ; ALFA_RBOUND = +10.001
          END IF

            CALL MR_GEN_INI_DZB_PROFILE_COEFFIS( NI , DZB_BK_L , DZB_BK_R , RTH_BK_L , RTH_BK_R , RTH , WTH ,   &
            & ALFA_LBOUND , ALFA_RBOUND , ALFA , BETA , GAMA , ERROR , ERRMSG )
            IF( ERROR < 0 ) THEN
              DEALLOCATE( ALFA , BETA , GAMA )
              RETURN
            END IF

            DYJ = 1.0_PARD_KIND / NJ

            DO J = 1 , NJ
             !DIR$ VECTOR ALIGNED
              DO I = 1 , NI
                ZB( I , J ) = ( GAMA( I ) - BETA( I ) * EXP( -ALFA( I ) * ( 1.0 - (J-0.5)*DYJ ) ) ) - HTH
              END DO
            END DO

          DEALLOCATE( ALFA_LBOUND , ALFA_RBOUND )

        DEALLOCATE( ALFA , BETA , GAMA )

      DEALLOCATE( WTH , RTH , RTH_BK_L , RTH_BK_R )

    DEALLOCATE( DZB_BK_L , DZB_BK_R )

  END SUBROUTINE MR_GEN_INI_ZB

!***********************************************************************************************************************************
! UNIT:
!
!  (SUBROUTINE)
!
! PURPOSE:
!
!   TO
!
! DEFINITION OF VARIABLES:
!
!
!
! RECORD OF REVISIONS:
!
!      DATE       |    PROGRAMMER    |    DESCRIPTION OF CHANGE
!      ====       |    ==========    |    =====================
!   20XX-XX-XX    |     DR. HYDE     |    ORIGINAL CODE.
!
!***********************************************************************************************************************************
  SUBROUTINE MR_GEN_INI_DZB_BK( XI0 , COEFFI_XXIM , COEFFI_DZB_BK_MIN_N_MAX_1 , COEFFI_DZB_BK_MIN_N_MAX_2 ,   &
  & NBENDS , NI , DZB_BK_L , DZB_BK_R , ERROR , ERRMSG )

    IMPLICIT NONE

    REAL   (PARD_KIND) , INTENT(IN ) :: XI0
    REAL   (PARD_KIND) , INTENT(IN ) :: COEFFI_XXIM
    REAL   (PARD_KIND) , INTENT(IN ) :: COEFFI_DZB_BK_MIN_N_MAX_1 , COEFFI_DZB_BK_MIN_N_MAX_2

    INTEGER            , INTENT(IN ) :: NBENDS

    INTEGER(IJID_KIND) , INTENT(IN ) :: NI

    REAL   (PARD_KIND) , INTENT(OUT) , DIMENSION(1:NI1(NI,PARD_KIND)) :: DZB_BK_L , DZB_BK_R

    INTEGER            , INTENT(OUT) :: ERROR
    CHARACTER(   *   ) , INTENT(OUT) :: ERRMSG

    REAL   (PARD_KIND) , DIMENSION(1:NI1(NI,PARD_KIND)) :: XXI , XXXI

    REAL   (PARD_KIND) :: DUMMY_DZB_BK_1 , DUMMY_DZB_BK_2

    REAL   (PARD_KIND) :: DXI

    CHARACTER( 2**03 ) :: I_CHAR
    INTEGER(IJID_KIND) :: I

    DXI = 0.5_PARD_KIND * NBENDS / NI

   !DIR$ VECTOR ALIGNED
    DO I = 1 , NI
      XXI( I ) = (I-0.5)*DXI - XI0
    END DO

   !DIR$ NOVECTOR
    DO I = 1 , NI
      CALL MR_NEWTON_SOLVE_FEQ_XXXI( COEFFI_XXIM , XXI( I ) , XXXI( I ) , XXI( I ) , ERROR , ERRMSG )
      IF( ERROR < 0 ) THEN
        WRITE( I_CHAR , '(I<LEN(I_CHAR)>)' ) I
        ERRMSG = TRIM(ERRMSG)//" for cross section no. "//TRIM(ADJUSTL(I_CHAR))
        RETURN
      END IF
    END DO

   !DIR$ VECTOR ALIGNED
    DO I = 1 , NI

      DUMMY_DZB_BK_1 = COEFFI_DZB_BK_MIN_N_MAX_1 - COEFFI_DZB_BK_MIN_N_MAX_1 * COS( 4.0 * PI * XXXI( I ) )
      DUMMY_DZB_BK_2 = COEFFI_DZB_BK_MIN_N_MAX_2 * SIN( 2.0 * PI * XXXI( I ) )

      DZB_BK_L( I ) = DUMMY_DZB_BK_1 - DUMMY_DZB_BK_2
      DZB_BK_R( I ) = DUMMY_DZB_BK_1 + DUMMY_DZB_BK_2

    END DO

  END SUBROUTINE MR_GEN_INI_DZB_BK

!***********************************************************************************************************************************
! UNIT:
!
!  (SUBROUTINE)
!
! PURPOSE:
!
!   TO
!
! DEFINITION OF VARIABLES:
!
!
!
! RECORD OF REVISIONS:
!
!      DATE       |    PROGRAMMER    |    DESCRIPTION OF CHANGE
!      ====       |    ==========    |    =====================
!   20XX-XX-XX    |     DR. HYDE     |    ORIGINAL CODE.
!
!***********************************************************************************************************************************
  SUBROUTINE MR_NEWTON_SOLVE_FEQ_XXXI( COEFFI_XXIM , XXI , XXXI , XXXI_NEAR , ERROR , ERRMSG )

    IMPLICIT NONE

    REAL   (PARD_KIND) , INTENT(IN ) :: COEFFI_XXIM

    REAL   (PARD_KIND) , INTENT(IN ) :: XXI

    REAL   (PARD_KIND) , INTENT(OUT) :: XXXI

    REAL   (PARD_KIND) , INTENT(IN ) :: XXXI_NEAR

    REAL   (PARD_KIND) , PARAMETER   :: EPS_XXXI = 10.0*EPSILON(EPS_XXXI)
    REAL   (CARD_KIND) , PARAMETER   :: EPS_FEQ_XXXI = 10.0*EPSILON(EPS_FEQ_XXXI)
    REAL   (CARD_KIND) , PARAMETER   :: EPS_FFEQ_XXXI = 10.0*EPSILON(EPS_FFEQ_XXXI)

    INTEGER            , PARAMETER   :: N_ITERS = 10000

    INTEGER            , INTENT(OUT) :: ERROR
    CHARACTER(   *   ) , INTENT(OUT) :: ERRMSG

    REAL   (PARD_KIND) :: XXXI0 , XXXI1
    REAL   (CARD_KIND) :: FEQ_XXXI0 , FEQ_XXXI1
    REAL   (CARD_KIND) :: FFEQ_XXXI0

    INTEGER            :: I_ITER

    ERRMSG = ""

    XXXI0 = XXXI_NEAR
    FEQ_XXXI0 = MR_FUNC_FEQ_XXXI( COEFFI_XXIM , XXI , XXXI0 )

    DO I_ITER = 1 , N_ITERS

      FFEQ_XXXI0 = MR_FUNC_FFEQ_XXXI( COEFFI_XXIM , XXI , XXXI0 )
      IF( FFEQ_XXXI0 <= EPS_FFEQ_XXXI ) THEN
        XXXI = HUGE(XXXI)
        ERROR = ERROR_NEWTON_SOLVE_ZERO_DERIVATIVE
        ERRMSG = "Zero deivative when using newton method to solve xxxi"
        RETURN
      END IF

      XXXI1 = XXXI0 - FEQ_XXXI0 / FFEQ_XXXI0

      FEQ_XXXI1 = MR_FUNC_FEQ_XXXI( COEFFI_XXIM , XXI , XXXI1 )
      IF( ABS( XXXI1 - XXXI0 ) <= EPS_XXXI .OR. ABS( FEQ_XXXI1 ) <= EPS_FEQ_XXXI ) THEN
        XXXI = XXXI1
        ERROR = 0
        RETURN
      ELSE
        XXXI0 = XXXI1
        FEQ_XXXI0 = FEQ_XXXI1
      END IF

    END DO

    XXXI = HUGE(XXXI)
    ERROR = ERROR_NEWTON_SOLVE_MAX_NUMBER_OF_ITERATION
    ERRMSG = "Max number of iterations when using newton method to solve xxxi"

  END SUBROUTINE MR_NEWTON_SOLVE_FEQ_XXXI

!***********************************************************************************************************************************
! UNIT:
!
!  (SUBROUTINE)
!
! PURPOSE:
!
!   TO
!
! DEFINITION OF VARIABLES:
!
!
!
! RECORD OF REVISIONS:
!
!      DATE       |    PROGRAMMER    |    DESCRIPTION OF CHANGE
!      ====       |    ==========    |    =====================
!   20XX-XX-XX    |     DR. HYDE     |    ORIGINAL CODE.
!
!***********************************************************************************************************************************
  FUNCTION MR_FUNC_FEQ_XXXI( COEFFI_XXIM , XXI , XXXI ) RESULT( FEQ_XXXI )

   !DIR$ ATTRIBUTES VECTOR : UNIFORM( COEFFI_XXIM ) :: MR_FUNC_FEQ_XXXI

    IMPLICIT NONE

    REAL   (PARD_KIND) , INTENT(IN ) :: COEFFI_XXIM

    REAL   (PARD_KIND) , INTENT(IN ) :: XXI
    REAL   (PARD_KIND) , INTENT(IN ) :: XXXI

    REAL   (CARD_KIND) :: FEQ_XXXI

    FEQ_XXXI = (XXXI-XXI) - COEFFI_XXIM * ( 1.0 - COS( 2.0 * PI * (XXXI+XXI) ) )

  END FUNCTION MR_FUNC_FEQ_XXXI

!***********************************************************************************************************************************
! UNIT:
!
!  (SUBROUTINE)
!
! PURPOSE:
!
!   TO
!
! DEFINITION OF VARIABLES:
!
!
!
! RECORD OF REVISIONS:
!
!      DATE       |    PROGRAMMER    |    DESCRIPTION OF CHANGE
!      ====       |    ==========    |    =====================
!   20XX-XX-XX    |     DR. HYDE     |    ORIGINAL CODE.
!
!***********************************************************************************************************************************
  FUNCTION MR_FUNC_FFEQ_XXXI( COEFFI_XXIM , XXI , XXXI ) RESULT( FFEQ_XXXI )

   !DIR$ ATTRIBUTES VECTOR : UNIFORM( COEFFI_XXIM ) :: MR_FUNC_FFEQ_XXXI

    IMPLICIT NONE

    REAL   (PARD_KIND) , INTENT(IN ) :: COEFFI_XXIM

    REAL   (PARD_KIND) , INTENT(IN ) :: XXI
    REAL   (PARD_KIND) , INTENT(IN ) :: XXXI

    REAL   (CARD_KIND) :: FFEQ_XXXI

    FFEQ_XXXI = 1.0 - 2.0 * PI * COEFFI_XXIM * SIN( 2.0 * PI * (XXXI+XXI) )

  END FUNCTION MR_FUNC_FFEQ_XXXI

!***********************************************************************************************************************************
! UNIT:
!
!  (SUBROUTINE)
!
! PURPOSE:
!
!   TO
!
! DEFINITION OF VARIABLES:
!
!
!
! RECORD OF REVISIONS:
!
!      DATE       |    PROGRAMMER    |    DESCRIPTION OF CHANGE
!      ====       |    ==========    |    =====================
!   20XX-XX-XX    |     DR. HYDE     |    ORIGINAL CODE.
!
!***********************************************************************************************************************************
  SUBROUTINE MR_GEN_INI_DZB_PROFILE_COEFFIS( NI , DZB_BK_L , DZB_BK_R , RTH_BK_L , RTH_BK_R , RTH , WTH ,   &
  & ALFA_LBOUND , ALFA_RBOUND , ALFA , BETA , GAMA , ERROR , ERRMSG )

    IMPLICIT NONE

    INTEGER(IJID_KIND) , INTENT(IN ) :: NI

    REAL   (PARD_KIND) , INTENT(IN ) , DIMENSION(1:NI1(NI,PARD_KIND)) :: DZB_BK_L , DZB_BK_R
    REAL   (GJRD_KIND) , INTENT(IN ) , DIMENSION(1:NI1(NI,GJRD_KIND)) :: RTH_BK_L , RTH_BK_R , RTH , WTH

    REAL   (CARD_KIND) , INTENT(IN ) , DIMENSION(1:NI1(NI,CARD_KIND)) :: ALFA_LBOUND , ALFA_RBOUND

    REAL   (CARD_KIND) , INTENT(OUT) , DIMENSION(1:NI1(NI,CARD_KIND)) :: ALFA , BETA , GAMA

    INTEGER            , INTENT(OUT) :: ERROR
    CHARACTER(   *   ) , INTENT(OUT) :: ERRMSG

    CHARACTER( 2**03 ) :: I_CHAR
    INTEGER(IJID_KIND) :: I

   !DIR$ NOVECTOR
    DO I = 1 , NI
      CALL MR_BISECT_SOLVE_FEQ_ALFA( DZB_BK_L( I ) , DZB_BK_R( I ) , RTH_BK_L( I ) , RTH_BK_R( I ) , RTH( I ) , WTH( I ) ,   &
      & ALFA( I ) , ALFA_LBOUND( I ) , ALFA_RBOUND( I ) , ERROR , ERRMSG )
      IF( ERROR < 0 ) THEN
        WRITE( I_CHAR , '(I<LEN(I_CHAR)>)' ) I
        ERRMSG = TRIM(ERRMSG)//" for cross section no. "//TRIM(ADJUSTL(I_CHAR))
        RETURN
      END IF
    END DO

   !DIR$ VECTOR ALIGNED
    DO I = 1 , NI
      CALL MR_CALC_BETA_N_GAMA( DZB_BK_L( I ) , DZB_BK_R( I ) , ALFA( I ) , BETA( I ) , GAMA( I ) )
    END DO

  END SUBROUTINE MR_GEN_INI_DZB_PROFILE_COEFFIS

!***********************************************************************************************************************************
! UNIT:
!
!  (SUBROUTINE)
!
! PURPOSE:
!
!   TO
!
! DEFINITION OF VARIABLES:
!
!
!
! RECORD OF REVISIONS:
!
!      DATE       |    PROGRAMMER    |    DESCRIPTION OF CHANGE
!      ====       |    ==========    |    =====================
!   20XX-XX-XX    |     DR. HYDE     |    ORIGINAL CODE.
!
!***********************************************************************************************************************************
  SUBROUTINE MR_BISECT_SOLVE_FEQ_ALFA( DZB_BK_L , DZB_BK_R , RTH_BK_L , RTH_BK_R , RTH , WTH ,   &
  & ALFA , ALFA_LBOUND , ALFA_RBOUND , ERROR , ERRMSG )

    IMPLICIT NONE

    REAL   (PARD_KIND) , INTENT(IN ) :: DZB_BK_L , DZB_BK_R
    REAL   (GJRD_KIND) , INTENT(IN ) :: RTH_BK_L , RTH_BK_R , RTH , WTH

    REAL   (CARD_KIND) , INTENT(OUT) :: ALFA

    REAL   (CARD_KIND) , INTENT(IN ) :: ALFA_LBOUND , ALFA_RBOUND

    REAL   (CARD_KIND) , PARAMETER   :: EPS_ALFA = 10.0*EPSILON(EPS_ALFA)
    REAL   (CARD_KIND) , PARAMETER   :: EPS_FEQ_ALFA = 10.0*EPSILON(EPS_FEQ_ALFA)

    INTEGER            , INTENT(OUT) :: ERROR
    CHARACTER(   *   ) , INTENT(OUT) :: ERRMSG

    REAL   (CARD_KIND) :: ALFA_LAPPRX , ALFA_RAPPRX , ALFA_MIDDLE
    REAL   (CARD_KIND) :: FEQ_ALFA_LAPPRX , FEQ_ALFA_RAPPRX , FEQ_ALFA_MIDDLE

    ERRMSG = ""

    ALFA_LAPPRX = ALFA_LBOUND
    FEQ_ALFA_LAPPRX = MR_FUNC_FEQ_ALFA( DZB_BK_L , DZB_BK_R , RTH_BK_L , RTH_BK_R , RTH , WTH , ALFA_LAPPRX )
    IF( ABS( FEQ_ALFA_LAPPRX ) <= EPS_FEQ_ALFA ) THEN
      ALFA = ALFA_LAPPRX
      ERROR = 0
      RETURN
    END IF
    ALFA_RAPPRX = ALFA_RBOUND
    FEQ_ALFA_RAPPRX = MR_FUNC_FEQ_ALFA( DZB_BK_L , DZB_BK_R , RTH_BK_L , RTH_BK_R , RTH , WTH , ALFA_RAPPRX )
    IF( ABS( FEQ_ALFA_RAPPRX ) <= EPS_FEQ_ALFA ) THEN
      ALFA = ALFA_RAPPRX
      ERROR = 0
      RETURN
    END IF

    IF( FEQ_ALFA_LAPPRX * FEQ_ALFA_RAPPRX > 0 ) THEN
      ALFA = HUGE(ALFA)
      ERROR = ERROR_BISECT_SOLVE_NO_UNIQUE_ROOT_IN_REGION
      ERRMSG = "No unique root in given region when using bisection method to solve alpha"
      RETURN
    END IF

    DO WHILE( 0.5 * ABS( ALFA_RAPPRX - ALFA_LAPPRX ) > EPS_ALFA )

      ALFA_MIDDLE = 0.5 * ( ALFA_LAPPRX + ALFA_RAPPRX )
      FEQ_ALFA_MIDDLE = MR_FUNC_FEQ_ALFA( DZB_BK_L , DZB_BK_R , RTH_BK_L , RTH_BK_R , RTH , WTH , ALFA_MIDDLE )

      IF( ABS( FEQ_ALFA_MIDDLE ) <= EPS_FEQ_ALFA ) THEN
        EXIT
      ELSE IF( FEQ_ALFA_LAPPRX * FEQ_ALFA_MIDDLE < 0 ) THEN
        ALFA_RAPPRX = ALFA_MIDDLE
        FEQ_ALFA_RAPPRX = FEQ_ALFA_MIDDLE
      ELSE
        ALFA_LAPPRX = ALFA_MIDDLE
        FEQ_ALFA_LAPPRX = FEQ_ALFA_MIDDLE
      END IF

    END DO

    ALFA = 0.5 * ( ALFA_LAPPRX + ALFA_RAPPRX )
    ERROR = 0

  END SUBROUTINE MR_BISECT_SOLVE_FEQ_ALFA

!***********************************************************************************************************************************
! UNIT:
!
!  (SUBROUTINE)
!
! PURPOSE:
!
!   TO
!
! DEFINITION OF VARIABLES:
!
!
!
! RECORD OF REVISIONS:
!
!      DATE       |    PROGRAMMER    |    DESCRIPTION OF CHANGE
!      ====       |    ==========    |    =====================
!   20XX-XX-XX    |     DR. HYDE     |    ORIGINAL CODE.
!
!***********************************************************************************************************************************
  FUNCTION MR_FUNC_FEQ_ALFA( DZB_BK_L , DZB_BK_R , RTH_BK_L , RTH_BK_R , RTH , WTH , ALFA ) RESULT( FEQ_ALFA )

   !DIR$ ATTRIBUTES VECTOR :: MR_FUNC_FEQ_ALFA

    IMPLICIT NONE

    REAL   (PARD_KIND) , INTENT(IN ) :: DZB_BK_L , DZB_BK_R
    REAL   (GJRD_KIND) , INTENT(IN ) :: RTH_BK_L , RTH_BK_R , RTH , WTH

    REAL   (CARD_KIND) , INTENT(IN ) :: ALFA

    REAL   (CARD_KIND) :: BETA , GAMA

    REAL   (CARD_KIND) :: FEQ_ALFA

    CALL MR_CALC_BETA_N_GAMA( DZB_BK_L , DZB_BK_R , ALFA , BETA , GAMA )

    FEQ_ALFA =   &
    & BETA * ( EXP(-ALFA) * ( ALFA * ( RTH_BK_R / WTH ) - 1.0 ) - ALFA * ( RTH_BK_L / WTH ) + 1.0 ) +   &
    & GAMA * ALFA * ALFA * ( RTH / WTH )

  END FUNCTION MR_FUNC_FEQ_ALFA

!***********************************************************************************************************************************
! UNIT:
!
!  (SUBROUTINE)
!
! PURPOSE:
!
!   TO
!
! DEFINITION OF VARIABLES:
!
!
!
! RECORD OF REVISIONS:
!
!      DATE       |    PROGRAMMER    |    DESCRIPTION OF CHANGE
!      ====       |    ==========    |    =====================
!   20XX-XX-XX    |     DR. HYDE     |    ORIGINAL CODE.
!
!***********************************************************************************************************************************
  SUBROUTINE MR_CALC_BETA_N_GAMA( DZB_BK_L , DZB_BK_R , ALFA , BETA , GAMA )

   !DIR$ ATTRIBUTES VECTOR :: MR_CALC_BETA_N_GAMA

    IMPLICIT NONE

    REAL   (PARD_KIND) , INTENT(IN ) :: DZB_BK_L , DZB_BK_R

    REAL   (CARD_KIND) , INTENT(IN ) :: ALFA

    REAL   (CARD_KIND) , INTENT(OUT) :: BETA , GAMA

    BETA = ( DZB_BK_L - DZB_BK_R ) / ( EXP(-ALFA) - 1.0 )
    GAMA =   DZB_BK_L + BETA

  END SUBROUTINE MR_CALC_BETA_N_GAMA

  END MODULE MR_MOD_GEN_INI_ZB


text:
! ****************************************************
!   output.f for sden2dRz
!  09 Feb. 2017    written by D.KAWATA
! ****************************************************

subroutine output(step,nxgrid,nygrid,xrange,yrange,tu,flagc)
      use gcdp_const
      use gcdp_system
      use mesh

      implicit none
      include 'mpif.h'

      integer,intent(in) :: step,nxgrid,nygrid,flagc
      double precision,intent(in) :: xrange(0:1),yrange(0:1),tu

      integer i,j,nval,ic,nmp,nmop,ip,npj
      character fileo*60
      double precision tym,tzmi,angc
! for work
      integer ierr
      integer istatus(MPI_STATUS_SIZE)
      integer,allocatable :: npjr(:),ireqs(:),ireqr(:)
      double precision,allocatable :: tdvs(:),tdvr(:)

      if(myrank.eq.0) then
        if(flagc.eq.0) then
          write(fileo,'(a15,i6.6)') 'output/gden2dRz',step
        else if(flagc.eq.1) then
          write(fileo,'(a15,i6.6)') 'output/dden2dRz',step
        else
          write(fileo,'(a15,i6.6)') 'output/sden2dRz',step
        endif
        write(6,*) 'output file name=',fileo
        open(61,file=fileo,status='unknown',form='unformatted')
        write(61) nxgrid,nygrid
        write(61) xrange(0),xrange(1)
        write(61) yrange(0),yrange(1)
        write(61) tu*TMUGYR
! for ascii data
!        if(flagc.eq.0) then
!          write(fileo,'(a15,i6.6,a4)') 'output/gden2dRz',step,'.asc'
!        else if(flagc.eq.1) then
!          write(fileo,'(a15,i6.6,a4)') 'output/dden2dRz',step,'.asc'
!        else
!          write(fileo,'(a15,i6.6,a4)') 'output/sden2dRz',step,'.asc'
!        endif
!        write(6,*) 'output ascii file name=',fileo
!        open(60,file=fileo,status='unknown')
!        write(60,'(a1,2I10)') '#',nxgrid,nygrid
!        write(60,'(a1,2(1pE13.5))') '#',xrange(0),xrange(1)
!        write(60,'(a1,2(1pE13.5))') '#',yrange(0),yrange(1)
!        write(60,'(a1,1pE13.5)') '#',tu*TMUGYR

        do j=0,nygrid-1
          do i=0,nxgrid-1
            write(61) x_m(i,j),y_m(i,j),denxy_m(i,j),met_m(i,j)
          enddo
        enddo

!        do j=0,nygrid-1
!          do i=0,nxgrid-1
!            write(60,'(4(1pE13.5))') x_m(i,j),y_m(i,j) &
!             ,denxy_m(i,j)*MUSM/(1.0e6),met_m(i,j)
!          enddo
!        enddo
!        close(60)
        close(61)

      endif

end subroutine

max_stars_repo_name: kawatadaisuke/gcdp-ana

max_issues_repo_name: kawatadaisuke/gcdp-ana

max_forks_repo_name: kawatadaisuke/gcdp-ana


text:
      SUBROUTINE  BMATVEC01( M, N, P, IX, AX, V, Y )

      !!***********************************************************************
      !! Version "$Id: bmatvec.f 1 2017-06-10 18:05:20Z coats $"
      !! EDSS/Models-3 I/O API.
      !! Copyright (C) 1992-2002 MCNC and Carlie J. Coats, Jr.,
      !! (C) 2003-2010 Baron Advanced Meteorological Systems, and
      !! (C) 2014 UNC Institute for the Environment.
      !! Distributed under the GNU LESSER GENERAL PUBLIC LICENSE version 2.1
      !! See file "LGPL.txt" for conditions of use.
      !!.........................................................................
      !!  subroutine BMATVEC01 body starts at line  92:  non-layered single-indexed inputs and outputs
      !!  subroutine BMATVEC02      starts at line 113:  non-layered single-indexed inputs, double-indexed outputs
      !!  subroutine BMATVEC11      starts at line 163:  single-indexed inputs and outputs
      !!  subroutine BMATVEC12      starts at line 275:  single-indexed inputs, double-indexed outputs
      !!  subroutine BMATVEC21      starts at line 330:  double-indexed inputs, single-indexed output
      !!  subroutine BMATVEC22      starts at line 382:  double-indexed inputs and output
      !!  subroutine BMATVEC        starts at line 437:  fall-back for non-"USE M3UTILIO"
      !!
      !!  FUNCTION:  apply a 4-band sparse matrix to an array ("layered vector")
      !!             as used for bilinear interpolation of meteorology in LAYPOINT
      !!
      !!  NOTE:  Transposes from I/O subscript order V( M,P )
      !!         to computational subscript order    C( P,N )
      !!
      !!       For bilinear interpolation of gridded data having dimension NC,NR
      !!       to a list of locations having grid-normal coordinates <X(S),Y(S)>:
      !!       let C(S) and R(S) be INT( X(S) ) and INT( Y(S) ), P(S) and Q(S)
      !!       be AMOD( X(S), 1.0 ) and AMOD( Y(S), 1.0 ).  Then IX has the
      !!       single-indexing subscripts into the grid for the cell=corners
      !!       surrounding the <X(S),Y(S)>
      !!           IX(1,S) = C + NC * ( R - 1 ) + 1
      !!           IX(2,S) = C + NC * ( R - 1 )
      !!           IX(3,S) = C + NC *   R
      !!           IX(4,S) = C + NC *   R       + 1
      !!       and AX has the bilinear-interpolation coefficients:
      !!           AX(1,S) = ( 1.0 - P( S ) )*( 1.0 - Q( S ) )
      !!           AX(2,S) =         P( S )  *( 1.0 - Q( S ) )
      !!           AX(3,S) = ( 1.0 - P( S ) )*        Q( S )
      !!           AX(4,S) =         P( S )  *        Q( S )
      !!
      !!  SEE ALSO:
      !!       UNGRIDB() which produces such matrices
      !!       BILIN()   which performs combined interpolate-only,
      !!                 preserving the subscript-order.
      !!
      !!  PRECONDITIONS REQUIRED:
      !!        USE M3UTILIO for generic INTERFACE
      !!        Number of layers same for input and output.
      !!        Index and coefficients set up so that the equation for
      !!        the multiplication (at each row R and layer L) is
      !!
      !!       Y(L,R) = SUM_{ J=1...4 }  A( J,R ) * V( I( J,R ),L )
      !!
      !!  SUBROUTINES AND FUNCTIONS CALLED:  none
      !!
      !!  REVISION  HISTORY:
      !!       prototype 12/1995 by CJC
      !!       Version    9/2014 by CJC:  modifications for OpenMP parallel
      !!        Version  12/2014 by CJC for I/O API v3.2:  multiple versions
      !!        with  M3UTILIO generic interface BMATVEC()
      !!    Version  03/2016 by CJC:  Add BMATVEC() for backwards compatibility
      !!***********************************************************************

        IMPLICIT NONE

        !!...........   ARGUMENTS and their descriptions:

        INTEGER, INTENT(IN   ) :: M               ! length of input  vector
        INTEGER, INTENT(IN   ) :: N               ! length of output vector
        INTEGER, INTENT(IN   ) :: P               ! number of layers
        INTEGER, INTENT(IN   ) :: IX( 4,N )       ! index array
        REAL   , INTENT(IN   ) :: AX( 4,N )       ! 4-band coeff matrix
        REAL   , INTENT(IN   ) :: V( M )          ! non-layered input  vector
        REAL   , INTENT(  OUT) :: Y( N )          ! non-layered output vector

        !!...........   SCRATCH LOCAL VARIABLES and their descriptions:

        INTEGER         S, J1, J2, J3, J4

        !!***********************************************************************
        !!   begin body of subroutine  BMATVEC01

!$OMP   PARALLEL DO
!$OMP&    DEFAULT( NONE ),
!$OMP&     SHARED( N, P, IX, AX, V, Y ),
!$OMP&    PRIVATE( S, J1, J2, J3, J4 )

        DO  S = 1, N

            J1 = IX( 1,S )
            J2 = IX( 2,S )
            J3 = IX( 3,S )
            J4 = IX( 4,S )
            Y( S ) = AX( 1,S ) * V( J1 )  +
     &               AX( 2,S ) * V( J2 )  +
     &               AX( 3,S ) * V( J3 )  +
     &               AX( 4,S ) * V( J4 )

        END DO

        RETURN

      END SUBROUTINE BMATVEC01

!!-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

      SUBROUTINE  BMATVEC02( M, NC, NR, P, IX, AX, V, Y )

        IMPLICIT NONE

        !!...........   ARGUMENTS and their descriptions:

        INTEGER, INTENT(IN   ) :: M                 ! length of input  vector
        INTEGER, INTENT(IN   ) :: NC, NR            ! dims of output grid
        INTEGER, INTENT(IN   ) :: P                 ! number of layers
        INTEGER, INTENT(IN   ) :: IX( 4,NC*NR )     ! index array
        REAL   , INTENT(IN   ) :: AX( 4,NC*NR )     ! 4-band coeff matrix
        REAL   , INTENT(IN   ) :: V( M  )           ! non-layered input  vector
        REAL   , INTENT(  OUT) :: Y( NC,NR )        ! non-layered output vector

        !!...........   SCRATCH LOCAL VARIABLES and their descriptions:

        INTEGER         C, R, S, J1, J2, J3, J4

        !!***********************************************************************
        !!   begin body of subroutine  BMATVEC01

!$OMP   PARALLEL DO
!$OMP&    DEFAULT( NONE ),
!$OMP&     SHARED( NC, NR, P, IX, AX, V, Y ),
!$OMP&    PRIVATE( C, R, S, J1, J2, J3, J4 )

        DO  R = 1, NR
        DO  C = 1, NC

            S = ( R  - 1 ) * NC  +  C
            J1 = IX( 1,S )
            J2 = IX( 2,S )
            J3 = IX( 3,S )
            J4 = IX( 4,S )
            Y( C,R ) = AX( 1,S ) * V( J1 )  +
     &                 AX( 2,S ) * V( J2 )  +
     &                 AX( 3,S ) * V( J3 )  +
     &                 AX( 4,S ) * V( J4 )

        END DO
        END DO

        RETURN

      END SUBROUTINE BMATVEC02

      SUBROUTINE  BMATVEC11( M, N, P, IX, AX, V, Y )

      !!***********************************************************************
      !! Version "$Id: bmatvec.f 1 2017-06-10 18:05:20Z coats $"
      !! EDSS/Models-3 I/O API.
      !! Copyright (C) 1992-2002 MCNC and Carlie J. Coats, Jr.,
      !! (C) 2003-2010 Baron Advanced Meteorological Systems, and
      !! (C) 2014 UNC Institute for the Environment.
      !! Distributed under the GNU LESSER GENERAL PUBLIC LICENSE version 2.1
      !! See file "LGPL.txt" for conditions of use.
      !!.........................................................................
      !!  subroutine BMATVEC11 body starts at line  90:  single-indexed inputs and outputs
      !!  subroutine BMATVEC12 body starts at line 114:  single-indexed inputs, double-indexed outputs
      !!  subroutine BMATVEC21 body starts at line 169:  double-indexed inputs, single-indexed output
      !!  subroutine BMATVEC22 body starts at line 221:  double-indexed inputs and output
      !!  subroutine BMATVEC   body starts at line 276:  fall-back for non-"USE M3UTILIO"
      !!
      !!  FUNCTION:  apply a 4-band sparse matrix to an array ("layered vector")
      !!             as used for bilinear interpolation of meteorology in LAYPOINT
      !!
      !!  NOTE:  Transposes from I/O subscript order V( M,P )
      !!         to computational subscript order    C( P,N )
      !!
      !!       For bilinear interpolation of gridded data having dimension NC,NR
      !!       to a list of locations having grid-normal coordinates <X(S),Y(S)>:
      !!       let C(S) and R(S) be INT( X(S) ) and INT( Y(S) ), P(S) and Q(S)
      !!       be AMOD( X(S), 1.0 ) and AMOD( Y(S), 1.0 ).  Then IX has the
      !!       single-indexing subscripts into the grid for the cell=corners
      !!       surrounding the <X(S),Y(S)>
      !!           IX(1,S) = C + NC * ( R - 1 ) + 1
      !!           IX(2,S) = C + NC * ( R - 1 )
      !!           IX(3,S) = C + NC *   R
      !!           IX(4,S) = C + NC *   R       + 1
      !!       and AX has the bilinear-interpolation coefficients:
      !!           AX(1,S) = ( 1.0 - P( S ) )*( 1.0 - Q( S ) )
      !!           AX(2,S) =         P( S )  *( 1.0 - Q( S ) )
      !!           AX(3,S) = ( 1.0 - P( S ) )*        Q( S )
      !!           AX(4,S) =         P( S )  *        Q( S )
      !!
      !!  SEE ALSO:
      !!       UNGRIDB() which produces such matrices
      !!       BILIN()   which performs combined interpolate-only,
      !!                 preserving the subscript-order.
      !!
      !!  PRECONDITIONS REQUIRED:
      !!        USE M3UTILIO for generic INTERFACE
      !!        Number of layers same for input and output.
      !!        Index and coefficients set up so that the equation for
      !!        the multiplication (at each row R and layer L) is
      !!
      !!       Y(L,R) = SUM_{ J=1...4 }  A( J,R ) * V( I( J,R ),L )
      !!
      !!  SUBROUTINES AND FUNCTIONS CALLED:  none
      !!
      !!  REVISION  HISTORY:
      !!       prototype 12/1995 by CJC
      !!       Version    9/2014 by CJC:  modifications for OpenMP parallel
      !!        Version  12/2014 by CJC for I/O API v3.2:  multiple versions
      !!        with  M3UTILIO generic interface BMATVEC()
      !!    Version  03/2016 by CJC:  Add BMATVEC() for backwards compatibility
      !!***********************************************************************

        IMPLICIT NONE

        !!...........   ARGUMENTS and their descriptions:

        INTEGER, INTENT(IN   ) :: M               ! length of input  vector
        INTEGER, INTENT(IN   ) :: N               ! length of output vector
        INTEGER, INTENT(IN   ) :: P               ! number of layers
        INTEGER, INTENT(IN   ) :: IX( 4,N )       ! index array
        REAL   , INTENT(IN   ) :: AX( 4,N )       ! 4-band coeff matrix
        REAL   , INTENT(IN   ) :: V( M,P )        ! P-layered input  vector
        REAL   , INTENT(  OUT) :: Y( P,N )        ! P-layered output vector

        !!...........   SCRATCH LOCAL VARIABLES and their descriptions:

        INTEGER         S, L, J1, J2, J3, J4

        !!***********************************************************************
        !!   begin body of subroutine  BMATVEC11

!$OMP   PARALLEL DO
!$OMP&    DEFAULT( NONE ),
!$OMP&     SHARED( N, P, IX, AX, V, Y ),
!$OMP&    PRIVATE( S, L, J1, J2, J3, J4 )

        DO  S = 1, N

            J1 = IX( 1,S )
            J2 = IX( 2,S )
            J3 = IX( 3,S )
            J4 = IX( 4,S )

            DO  L = 1, P
                Y( L,S ) = AX( 1,S ) * V( J1,L )  +
     &                     AX( 2,S ) * V( J2,L )  +
     &                     AX( 3,S ) * V( J3,L )  +
     &                     AX( 4,S ) * V( J4,L )
            END DO

        END DO

        RETURN

      END SUBROUTINE BMATVEC11

!!-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

      SUBROUTINE  BMATVEC12( M, NC, NR, P, IX, AX, V, Y )

        IMPLICIT NONE

        !!...........   ARGUMENTS and their descriptions:

        INTEGER, INTENT(IN   ) :: M                 ! length of input  vector
        INTEGER, INTENT(IN   ) :: NC, NR            ! dims of output grid
        INTEGER, INTENT(IN   ) :: P                 ! number of layers
        INTEGER, INTENT(IN   ) :: IX( 4,NC*NR )     ! index array
        REAL   , INTENT(IN   ) :: AX( 4,NC*NR )     ! 4-band coeff matrix
        REAL   , INTENT(IN   ) :: V( M,P )          ! P-layered input  vector
        REAL   , INTENT(  OUT) :: Y( P,NC,NR )      ! P-layered output vector

        !!...........   SCRATCH LOCAL VARIABLES and their descriptions:

        INTEGER         C, R, S, L, J1, J2, J3, J4

        !!***********************************************************************
        !!   begin body of subroutine  BMATVEC11

!$OMP   PARALLEL DO
!$OMP&    DEFAULT( NONE ),
!$OMP&     SHARED( NC, NR, P, IX, AX, V, Y ),
!$OMP&    PRIVATE( C, R, S, L, J1, J2, J3, J4 )

        DO  R = 1, NR
        DO  C = 1, NC

            S = ( R  - 1 ) * NC  +  C
            J1 = IX( 1,S )
            J2 = IX( 2,S )
            J3 = IX( 3,S )
            J4 = IX( 4,S )

            DO  L = 1, P
                Y( L,C,R ) = AX( 1,S ) * V( J1,L )  +
     &                       AX( 2,S ) * V( J2,L )  +
     &                       AX( 3,S ) * V( J3,L )  +
     &                       AX( 4,S ) * V( J4,L )
            END DO

        END DO
        END DO

        RETURN

      END SUBROUTINE BMATVEC12

!!-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

      SUBROUTINE  BMATVEC21( MC, MR, N, P, IX, AX, V, Y )

        IMPLICIT NONE

        !!...........   ARGUMENTS and their descriptions:

        INTEGER, INTENT(IN   ) :: MC, MR            ! length of input  vector
        INTEGER, INTENT(IN   ) :: N                 ! dims of output grid
        INTEGER, INTENT(IN   ) :: P                 ! number of layers
        INTEGER, INTENT(IN   ) :: IX( 4,N )         ! index array
        REAL   , INTENT(IN   ) :: AX( 4,N )         ! 4-band coeff matrix
        REAL   , INTENT(IN   ) :: V( MC,MR,P )      ! P-layered input  vector
        REAL   , INTENT(  OUT) :: Y( P,N )          ! P-layered output vector

        !!...........   SCRATCH LOCAL VARIABLES and their descriptions:

        INTEGER         C, R, S, CC, RR, K, L

        !!***********************************************************************
        !!   begin body of subroutine  BMATVEC11

!$OMP   PARALLEL DO
!$OMP&    DEFAULT( NONE ),
!$OMP&     SHARED( N, MC, MR, P, IX, AX, V, Y ),
!$OMP&    PRIVATE( C, R, S, K, L, CC, RR )

        DO  S = 1, N

            K  = IX( 1,S )
            CC = MOD( K, MC )
            RR = 1 + K / MC

            DO  L = 1, P
                Y( L,S ) = AX( 1,S ) * V( CC  ,RR  ,L )  +
     &                     AX( 2,S ) * V( CC+1,RR  ,L )  +
     &                     AX( 3,S ) * V( CC  ,RR+1,L )  +
     &                     AX( 4,S ) * V( CC+1,RR+1,L )
            END DO

        END DO

        RETURN

      END SUBROUTINE BMATVEC21

!!-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

      SUBROUTINE  BMATVEC22( MC, MR, NC, NR, P, IX, AX, V, Y )

        IMPLICIT NONE

        !!...........   ARGUMENTS and their descriptions:

        INTEGER, INTENT(IN   ) :: MC, MR            ! length of input  vector
        INTEGER, INTENT(IN   ) :: NC, NR            ! dims of output grid
        INTEGER, INTENT(IN   ) :: P                 ! number of layers
        INTEGER, INTENT(IN   ) :: IX( 4,NC*NR )     ! index array
        REAL   , INTENT(IN   ) :: AX( 4,NC*NR )     ! 4-band coeff matrix
        REAL   , INTENT(IN   ) :: V( MC,MR,P )      ! P-layered input  vector
        REAL   , INTENT(  OUT) :: Y( P,NC,NR )      ! P-layered output vector

        !!...........   SCRATCH LOCAL VARIABLES and their descriptions:

        INTEGER         C, R, S, CC, RR, K, L

        !!***********************************************************************
        !!   begin body of subroutine  BMATVEC11

!$OMP   PARALLEL DO
!$OMP&    DEFAULT( NONE ),
!$OMP&     SHARED( NC, NR, MC, MR, P, IX, AX, V, Y ),
!$OMP&    PRIVATE( C, R, S, K, L, CC, RR )

        DO  R = 1, NR
        DO  C = 1, NC

            S  = ( R  - 1 ) * NC  +  C
            K  = IX( 1,S )
            CC = MOD( K, MC )
            RR = 1 + K / MC

            DO  L = 1, P
                Y( L,C,R ) = AX( 1,S ) * V( CC  ,RR  ,L )  +
     &                       AX( 2,S ) * V( CC+1,RR  ,L )  +
     &                       AX( 3,S ) * V( CC  ,RR+1,L )  +
     &                       AX( 4,S ) * V( CC+1,RR+1,L )
            END DO

        END DO
        END DO

        RETURN

      END SUBROUTINE BMATVEC22

!!-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

      SUBROUTINE  BMATVEC021( MC, MR, N, P, IX, AX, V, Y )

        IMPLICIT NONE

        !!...........   ARGUMENTS and their descriptions:

        INTEGER, INTENT(IN   ) :: MC, MR            ! length of input  vector
        INTEGER, INTENT(IN   ) :: N                 ! dims of output grid
        INTEGER, INTENT(IN   ) :: P                 ! number of layers
        INTEGER, INTENT(IN   ) :: IX( 4,N )         ! index array
        REAL   , INTENT(IN   ) :: AX( 4,N )         ! 4-band coeff matrix
        REAL   , INTENT(IN   ) :: V( MC,MR )        ! P-layered input  vector
        REAL   , INTENT(  OUT) :: Y( N )            ! P-layered output vector

        !!...........   SCRATCH LOCAL VARIABLES and their descriptions:

        INTEGER         C, R, S, CC, RR, K

        !!***********************************************************************
        !!   begin body of subroutine  BMATVEC11

!$OMP   PARALLEL DO
!$OMP&    DEFAULT( NONE ),
!$OMP&     SHARED( N, MC, MR, P, IX, AX, V, Y ),
!$OMP&    PRIVATE( C, R, S, K, CC, RR )

        DO  S = 1, N

            K  = IX( 1,S )
            CC = MOD( K, MC )
            RR = 1 + K / MC
            Y( S ) = AX( 1,S ) * V( CC  ,RR   )  +
     &               AX( 2,S ) * V( CC+1,RR   )  +
     &               AX( 3,S ) * V( CC  ,RR+1 )  +
     &               AX( 4,S ) * V( CC+1,RR+1 )

        END DO

        RETURN

      END SUBROUTINE BMATVEC021

!!-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

      SUBROUTINE  BMATVEC022( MC, MR, NC, NR, P, IX, AX, V, Y )

        IMPLICIT NONE

        !!...........   ARGUMENTS and their descriptions:

        INTEGER, INTENT(IN   ) :: MC, MR            ! length of input  vector
        INTEGER, INTENT(IN   ) :: NC, NR            ! dims of output grid
        INTEGER, INTENT(IN   ) :: P                 ! number of layers
        INTEGER, INTENT(IN   ) :: IX( 4,NC*NR )     ! index array
        REAL   , INTENT(IN   ) :: AX( 4,NC*NR )     ! 4-band coeff matrix
        REAL   , INTENT(IN   ) :: V( MC,MR )        ! P-layered input  vector
        REAL   , INTENT(  OUT) :: Y( NC,NR )        ! P-layered output vector

        !!...........   SCRATCH LOCAL VARIABLES and their descriptions:

        INTEGER         C, R, S, CC, RR, K

        !!***********************************************************************
        !!   begin body of subroutine  BMATVEC11

!$OMP   PARALLEL DO
!$OMP&    DEFAULT( NONE ),
!$OMP&     SHARED( NC, NR, MC, MR, P, IX, AX, V, Y ),
!$OMP&    PRIVATE( C, R, S, K, CC, RR )

        DO  R = 1, NR
        DO  C = 1, NC

            S  = ( R  - 1 ) * NC  +  C
            K  = IX( 1,S )
            CC = MOD( K, MC )
            RR = 1 + K / MC
            Y( C,R ) = AX( 1,S ) * V( CC  ,RR   )  +
     &                 AX( 2,S ) * V( CC+1,RR   )  +
     &                 AX( 3,S ) * V( CC  ,RR+1 )  +
     &                 AX( 4,S ) * V( CC+1,RR+1 )

        END DO
        END DO

        RETURN

      END SUBROUTINE BMATVEC022

!!-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

      SUBROUTINE  BMATVEC( M, N, P, IX, AX, V, Y )

      IMPLICIT NONE

C...........   ARGUMENTS and their descriptions:

        INTEGER, INTENT(IN   ) :: M               ! length of input  vector
        INTEGER, INTENT(IN   ) :: N               ! length of output vector
        INTEGER, INTENT(IN   ) :: P               ! number of layers
        INTEGER, INTENT(IN   ) :: IX( 4,N )       ! index array
        REAL   , INTENT(IN   ) :: AX( 4,N )       ! 4-band coeff matrix
        REAL   , INTENT(IN   ) :: V( M,P )        ! P-layered input  vector
        REAL   , INTENT(  OUT) :: Y( P,N )        ! P-layered output vector

        !!................................................................

        CALL BMATVEC11( M, N, P, IX, AX, V, Y )

        RETURN
      END SUBROUTINE  BMATVEC

ext: f

max_stars_repo_name: mdhillmancmcl/TheWorldAvatar-CMCL-Fork

max_issues_repo_name: mdhillmancmcl/TheWorldAvatar-CMCL-Fork

max_forks_repo_name: mdhillmancmcl/TheWorldAvatar-CMCL-Fork


text:
program ProgramName
    DO i=1,6<caret>
end program ProgramName

max_stars_repo_name: LChernigovskaya/fortran-plugin

max_issues_repo_name: LChernigovskaya/fortran-plugin

max_forks_repo_name: LChernigovskaya/fortran-plugin


text:
!
! SA_MR
!
! Module containing forward, tangent-linear and adjoint subroutines
! for specific amount to/from mixing ratio units conversion.
!
! CREATION HISTORY:
!       Written by:     Paul van Delst, 06-Aug-2009
!                       paul.vandelst@noaa.gov
!

MODULE SA_MR

  ! -----------------
  ! Environment setup
  ! -----------------
  ! Modules used
  USE Type_Kinds, ONLY: fp
  USE Message_Handler
  USE Profile_Utility_Parameters, ONLY: ZERO, ONE, &
                                        MAX_N_MOLECULAR_SPECIES, &
                                        MW_DRYAIR, MOLECULAR_WEIGHT, &
                                        G_TO_KG
  ! Disable all implicit typing
  IMPLICIT NONE

  ! ------------
  ! Visibilities
  ! ------------
  PRIVATE
  PUBLIC :: SA_to_MR, SA_to_MR_TL, SA_to_MR_AD
  PUBLIC :: MR_to_SA, MR_to_SA_TL, MR_to_SA_AD

  ! -----------------
  ! Module parameters
  ! -----------------
  CHARACTER(*), PARAMETER :: MODULE_VERSION_ID = &

CONTAINS

!------------------------------------------------------------------------------
!:sdoc+:
! NAME:
!       SA_to_MR
!
! PURPOSE:
!       Elemental subroutine to convert gas concentrations from specific amount
!       in g/kg to mass mixing ratio in g/kg.
!
! CALLING SEQUENCE:
!       CALL SA_to_MR( Specific_Amount          , &  ! Input
!                      Mixing_Ratio             , &  ! Output
!                      Water_Vapor = Water_Vapor  )  ! Optional input
!
! INPUTS:
!       Specific_Amount:   Specific amount of gas.
!                          Must be > or = 0.0
!                          UNITS:      g/kg
!                          TYPE:       REAL(fp)
!                          DIMENSION:  Scalar or any rank
!                          ATTRIBUTES: INTENT(IN)
!
! OUTPUTS:
!       Mixing_Ratio:      Mass mixing ratio of gas.
!                          Set to 0.0 if input < 0.0
!                          UNITS:      g/kg
!                          TYPE:       REAL(fp)
!                          DIMENSION:  Same as Specific_Amount argument
!                          ATTRIBUTES: INTENT(OUT)
!
! OPTIONAL INPUT:
!       Water_Vapor:       Water vapor mass mixing ratio. If this argument is
!                          not supplied, the input SPECIFIC_AMOUNT argument is
!                          assumed to be water vapor.
!                          Output set to 0.0 if < 0.0
!                          UNITS:      g/kg
!                          TYPE:       REAL(fp)
!                          DIMENSION:  Same as Specific_Amount argument
!                          ATTRIBUTES: OPTIONAL, INTENT(IN)
!
! PROCEDURE:
!       The specific amount is defined as the ratio of the mass of gas
!       to the total mass of air and water vapor in a given volume:
!
!                  Mg
!         SA = ---------   .......(1)
!               Md + Mw
!
!       where Mg = mass of gas
!             Mw = mass of water vapor
!             Md = mass of dry air (including Mg).
!
!       The gas mass mixing ratio is defined as the ratio of the mass
!       of the gas to the mass of dry air in a given volume:
!
!                Mg
!         MRg = ----    .......(2)
!                Md
!
!       Rearranging (1) and substituting in (2) gives,
!
!          1      Md + Mw
!         ---- = ---------
!          SA       Mg
!
!                 Md     Mw
!              = ---- + ----
!                 Mg     Mg
!
!                 Md     Mw   Md                                      Md
!              = ---- + ----.----    ( multiplying the second term by ---- )
!                 Mg     Md   Mg                                      Md
!
!                  1
!              = ----- ( 1 + MRw )
!                 MRg
!
!       therefore,
!
!         MRg = SA ( 1 + MRw )
!
!       for input units of g/g or kg/kg. For input units of g/kg then,
!
!         MRg = SA ( 1 + 0.001*MRw )   .......(3)
!
!       If the input specific amount is for water vapor (specific humidity)
!       then (3) becomes,
!
!         MRw = SA ( 1 + 0.001*MRw )
!
!       i.e.
!
!         MRw = SA + 0.001*MRw*SA
!
!       i.e.
!
!         MRw( 1 - 0.001*SA ) = SA
!
!       therefore,
!
!                      SA
!         MRw = ------------------   .......(4)
!                ( 1 - 0.001*SA )
!
!
! CREATION HISTORY:
!       Written by:     Paul van Delst, 05-Dec-2003
!                       paul.vandelst@noaa.gov
!:sdoc-:
!------------------------------------------------------------------------------

  ELEMENTAL SUBROUTINE SA_to_MR( &
    Specific_Amount, &  ! Input
    Mixing_Ratio   , &  ! Output
    Water_Vapor      )  ! Optional Input
    ! Arguments
    REAL(fp),           INTENT(IN)  :: Specific_Amount
    REAL(fp),           INTENT(OUT) :: Mixing_Ratio
    REAL(fp), OPTIONAL, INTENT(IN)  :: Water_Vapor

    ! Error checks
    ! ...Zero output for -ve input
    IF ( Specific_Amount < ZERO ) THEN
      Mixing_Ratio = ZERO
      RETURN
    ENDIF

    ! Calculate mixing ratio in g/kg
    IF ( PRESENT( Water_Vapor ) ) THEN
      ! ...Zero output for -ve input
      IF ( Water_Vapor < ZERO ) THEN
        Mixing_Ratio = ZERO
        RETURN
      ENDIF
      ! ...Gas other than water vapour
      Mixing_Ratio = Specific_Amount * ( ONE + ( G_TO_KG * Water_Vapor ) )
    ELSE
      ! ...Water vapour
      Mixing_Ratio = Specific_Amount / ( ONE - ( G_TO_KG * Specific_Amount ) )
    END IF

  END SUBROUTINE SA_to_MR

!------------------------------------------------------------------------------
!:sdoc+:
! NAME:
!       SA_to_MR_TL
!
! PURPOSE:
!       Tangent-linear form of elemental subroutine to convert gas
!       concentrations from specific amount in g/kg to mass mixing
!       ratio in g/kg.
!
! CALLING SEQUENCE:
!       CALL SA_to_MR_TL( Specific_Amount   , &  ! FWD Input
!                         Specific_Amount_TL, &  ! TL  Input
!                         Mixing_Ratio_TL   , &  ! TL  Output
!                         Water_Vapor       , &  ! Optional FWD Input
!                         Water_Vapor_TL      )  ! Optional TL  Input
!
! INPUTS:
!       Specific_Amount:     Specific amount of gas.
!                            Must be > or = 0.0
!                            UNITS:      g/kg
!                            TYPE:       REAL(fp)
!                            DIMENSION:  Scalar or any rank
!                            ATTRIBUTES: INTENT(IN)
!
!       Specific_Amount_TL:  Tangent-linear specific amount of gas.
!                            UNITS:      g/kg
!                            TYPE:       REAL(fp)
!                            DIMENSION:  Same as Specific_Amount argument
!                            ATTRIBUTES: INTENT(IN)
!
! OUTPUTS:
!       Mixing_Ratio_TL:     Tangent-linear mass mixing ratio of gas.
!                            Set to 0.0 if input Specific_Amount < 0.0
!                            UNITS:      g/kg
!                            TYPE:       REAL(fp)
!                            DIMENSION:  Same as Specific_Amount argument
!                            ATTRIBUTES: INTENT(OUT)
!
! OPTIONAL INPUT:
!       Water_Vapor:         Water vapor mass mixing ratio. If this argument is
!                            not supplied, the input SPECIFIC_AMOUNT argument is
!                            assumed to be water vapor.
!                            Output set to 0.0 if < 0.0
!                            UNITS:      g/kg
!                            TYPE:       REAL(fp)
!                            DIMENSION:  Same as Specific_Amount argument
!                            ATTRIBUTES: OPTIONAL, INTENT(IN)
!
!       Water_Vapor_TL:      Tangent-linear water vapor mass mixing ratio. This
!                            argument is ignored if the optional Water_Vapor
!                            argument is not supplied also.
!                            UNITS:      g/kg
!                            TYPE:       REAL(fp)
!                            DIMENSION:  Same as Specific_Amount argument
!                            ATTRIBUTES: OPTIONAL, INTENT(IN)
!
! CREATION HISTORY:
!       Written by:     Paul van Delst, 14-Aug-2009
!                       paul.vandelst@noaa.gov
!:sdoc-:
!------------------------------------------------------------------------------

  ELEMENTAL SUBROUTINE SA_to_MR_TL( &
    Specific_Amount   , &  ! FWD Input
    Specific_Amount_TL, &  ! TL  Input
    Mixing_Ratio_TL   , &  ! TL  Output
    Water_Vapor       , &  ! Optional FWD Input
    Water_Vapor_TL      )  ! Optional TL  Input
    ! Arguments
    REAL(fp),           INTENT(IN)  :: Specific_Amount
    REAL(fp),           INTENT(IN)  :: Specific_Amount_TL
    REAL(fp),           INTENT(OUT) :: Mixing_Ratio_TL
    REAL(fp), OPTIONAL, INTENT(IN)  :: Water_Vapor
    REAL(fp), OPTIONAL, INTENT(IN)  :: Water_Vapor_TL

    ! Error checks
    ! ...Zero output for -ve input
    IF ( Specific_Amount < ZERO ) THEN
      Mixing_Ratio_TL = ZERO
      RETURN
    ENDIF

    ! Tangent-linear conversion
    IF ( PRESENT(Water_Vapor) ) THEN
      ! ...Zero output for -ve input
      IF ( Water_Vapor < ZERO ) THEN
        Mixing_Ratio_TL = ZERO
        RETURN
      ENDIF
      ! ...Gas other than water vapour
      Mixing_Ratio_TL = ((ONE + (G_TO_KG * Water_Vapor)) * Specific_Amount_TL)
      IF ( PRESENT(Water_Vapor_TL) ) THEN
        Mixing_Ratio_TL = Mixing_Ratio_TL + (G_TO_KG * Specific_Amount * Water_Vapor_TL)
      END IF
    ELSE
      ! ...Water vapour
      Mixing_Ratio_TL = Specific_Amount_TL / (ONE - (G_TO_KG * Specific_Amount))**2
    END IF

  END SUBROUTINE SA_to_MR_TL

!------------------------------------------------------------------------------
!:sdoc+:
! NAME:
!       SA_to_MR_AD
!
! PURPOSE:
!       Adjoint form of elemental subroutine to convert gas concentrations
!       from specific amount in g/kg to mass mixing ratio in g/kg.
!
! CALLING SEQUENCE:
!       CALL SA_to_MR_AD( Specific_Amount                , &  ! FWD Input
!                         Mixing_Ratio_AD                , &  ! AD  Input
!                         Specific_Amount_AD             , &  ! AD  Output
!                         Water_Vapor    = Water_Vapor   , &  ! Optional FWD Input
!                         Water_Vapor_AD = Water_Vapor_AD  )  ! Optional AD  Output
!
! INPUTS:
!       Specific_Amount:     Specific amount of gas.
!                            Must be > or = 0.0
!                            UNITS:      g/kg
!                            TYPE:       REAL(fp)
!                            DIMENSION:  Scalar or any rank
!                            ATTRIBUTES: INTENT(IN)
!
!       Mixing_Ratio_AD:     Adjoint mass mixing ratio of gas.
!                            *** SET TO ZERO ON EXIT ***
!                            UNITS:      g/kg
!                            TYPE:       REAL(fp)
!                            DIMENSION:  Same as Specific_Amount argument
!                            ATTRIBUTES: INTENT(IN OUT)
!
! OUTPUTS:
!       Specific_Amount_AD:  Adjoint specific amount of gas.
!                            *** MUST HAVE VALUE ON ENTRY ***
!                            UNITS:      g/kg
!                            TYPE:       REAL(fp)
!                            DIMENSION:  Same as Specific_Amount argument
!                            ATTRIBUTES: INTENT(IN OUT)
!
!
! OPTIONAL INPUT:
!       Water_Vapor:         Water vapor mass mixing ratio. If this argument is
!                            not supplied, the input SPECIFIC_AMOUNT argument is
!                            assumed to be water vapor.
!                            Output set to 0.0 if < 0.0
!                            UNITS:      g/kg
!                            TYPE:       REAL(fp)
!                            DIMENSION:  Same as Specific_Amount argument
!                            ATTRIBUTES: OPTIONAL, INTENT(IN)
!
! OPTIONAL OUTPUT:
!       Water_Vapor_AD:      Adjoint water vapor mass mixing ratio. This
!                            argument is ignored if the optional Water_Vapor
!                            argument is not supplied also.
!                            *** MUST HAVE VALUE ON ENTRY ***
!                            UNITS:      g/kg
!                            TYPE:       REAL(fp)
!                            DIMENSION:  Same as Specific_Amount argument
!                            ATTRIBUTES: OPTIONAL, INTENT(IN OUT)
!
! CREATION HISTORY:
!       Written by:     Paul van Delst, 14-Aug-2009
!                       paul.vandelst@noaa.gov
!:sdoc-:
!------------------------------------------------------------------------------

  ELEMENTAL SUBROUTINE SA_to_MR_AD( &
    Specific_Amount   , &  ! FWD Input
    Mixing_Ratio_AD   , &  ! AD  Input
    Specific_Amount_AD, &  ! AD  Output
    Water_Vapor       , &  ! Optional FWD Input
    Water_Vapor_AD      )  ! Optional AD  Output
    ! Arguments
    REAL(fp),           INTENT(IN)     :: Specific_Amount
    REAL(fp),           INTENT(IN OUT) :: Mixing_Ratio_AD
    REAL(fp),           INTENT(IN OUT) :: Specific_Amount_AD
    REAL(fp), OPTIONAL, INTENT(IN)     :: Water_Vapor
    REAL(fp), OPTIONAL, INTENT(IN OUT) :: Water_Vapor_AD

    ! Error checks
    ! ...Zero output for -ve input
    IF ( Specific_Amount < ZERO ) THEN
      Mixing_Ratio_AD    = ZERO
      Specific_Amount_AD = ZERO
      RETURN
    ENDIF

    ! Adjoint conversion
    IF ( PRESENT(Water_Vapor) ) THEN
      ! ...Zero output for -ve input
      IF ( Water_Vapor < ZERO ) THEN
        Mixing_Ratio_AD    = ZERO
        Specific_Amount_AD = ZERO
        RETURN
      ENDIF
      ! ...Gas other than water vapour
      IF ( PRESENT(Water_Vapor_AD) ) THEN
        Water_Vapor_AD = Water_Vapor_AD + (G_TO_KG * Specific_Amount * Mixing_Ratio_AD)
      END IF
      Specific_Amount_AD = Specific_Amount_AD + ((ONE + (G_TO_KG * Water_Vapor)) * Mixing_Ratio_AD)
    ELSE
      ! ...Water vapour
      Specific_Amount_AD = Specific_Amount_AD + (Mixing_Ratio_AD / (ONE - (G_TO_KG * Specific_Amount))**2)
    END IF
    Mixing_Ratio_AD = ZERO

  END SUBROUTINE SA_to_MR_AD

!------------------------------------------------------------------------------
!:sdoc+:
! NAME:
!       MR_to_SA
!
! PURPOSE:
!       Elemental subroutine to convert gas concentrations from mass mixing
!       ratio in g/kg to specific amount in g/kg.
!
! CALLING SEQUENCE:
!       CALL MR_to_SA( Mixing_Ratio             , &  ! Input
!                      Specific_Amount          , &  ! Output
!                      Water_Vapor = Water_Vapor  )  ! Optional input
!
! INPUTS:
!       Mixing_Ratio:      Mass mixing ratio of gas.
!                          Must be > or = 0.0
!                          UNITS:      g/kg
!                          TYPE:       REAL(fp)
!                          DIMENSION:  Scalar or any rank
!                          ATTRIBUTES: INTENT(IN)
!
! OUTPUTS:
!       Specific_Amount:   Specific amount of gas.
!                          Set to 0.0 if input < 0.0
!                          UNITS:      g/kg
!                          TYPE:       REAL(fp)
!                          DIMENSION:  Same as Mixing_Ratio argument
!                          ATTRIBUTES: INTENT(OUT)
!
! OPTIONAL INPUT:
!       Water_Vapor:       Water vapor specific humidity. If this argument is
!                          not supplied, the input MIXING_RATIO argument is
!                          assumed to be water vapor.
!                          Output set to 0.0 if < 0.0
!                          UNITS:      g/kg
!                          TYPE:       REAL(fp)
!                          DIMENSION:  Same as Mixing_Ratio argument
!                          ATTRIBUTES: OPTIONAL, INTENT(IN)
!
! PROCEDURE:
!       From the SA_to_MR conversion, we know that for input units of
!       g/kg that,
!
!         MRg = SAg ( 1 + 0.001*MRw )   .......(1)
!
!       and
!                      SAw
!         MRw = -------------------   .......(2)
!                ( 1 - 0.001*SAw )
!
!       where MRg = mass mixing ratio of gas
!             MRw = mass mixing ratio of water vapor
!             SAg = specific amount of gas
!             SAw = specific amount of water vapor (specific humidity)
!
!       Rearranging (1) and (2) gives,
!
!         SAg = MRg ( 1 - 0.001*SAw )
!
!       and for water vapor only,
!
!                       MRw
!         SAw = -------------------
!                ( 1 + 0.001*MRw )
!
! CREATION HISTORY:
!       Written by:     Paul van Delst, 05-Dec-2003
!                       paul.vandelst@noaa.gov
!:sdoc-:
!------------------------------------------------------------------------------

  ELEMENTAL SUBROUTINE MR_to_SA( &
    Mixing_Ratio   , &  ! Input
    Specific_Amount, &  ! Output
    Water_Vapor      )  ! Optional Input
    ! Arguments
    REAL(fp),           INTENT(IN)  :: Mixing_Ratio
    REAL(fp),           INTENT(OUT) :: Specific_Amount
    REAL(fp), OPTIONAL, INTENT(IN)  :: Water_Vapor

    ! Error checks
    ! ...Zero output for -ve input
    IF ( Mixing_Ratio < ZERO ) THEN
      Specific_Amount = ZERO
      RETURN
    ENDIF

    ! Calculate specific amount in g/kg
    IF ( PRESENT(Water_Vapor) ) THEN
      ! ...Zero output for -ve input
      IF ( Water_Vapor < ZERO ) THEN
        Specific_Amount = ZERO
        RETURN
      ENDIF
      ! ...Gas other than water vapour
      Specific_Amount = Mixing_Ratio * ( ONE - ( G_TO_KG * Water_Vapor ) )
    ELSE
      ! ...Water vapour
      Specific_Amount = Mixing_Ratio / ( ONE + ( G_TO_KG * Mixing_Ratio ) )
    END IF

  END SUBROUTINE MR_to_SA

!------------------------------------------------------------------------------
!:sdoc+:
! NAME:
!       MR_to_SA_TL
!
! PURPOSE:
!       Tangent-linear form of elemental subroutine to convert gas
!       concentrations from mass mixing ratio in g/kg to specific
!       amount in g/kg.
!
! CALLING SEQUENCE:
!       CALL MR_to_SA_TL( Mixing_Ratio                   , &  ! FWD Input
!                         Mixing_Ratio_TL                , &  ! TL  Input
!                         Specific_Amount_TL             , &  ! TL  Output
!                         Water_Vapor    = Water_Vapor   , &  ! Optional FWD input
!                         Water_Vapor_TL = Water_Vapor_TL  )  ! Optional TL input
!
! INPUTS:
!       Mixing_Ratio:        Mass mixing ratio of gas.
!                            Must be > or = 0.0
!                            UNITS:      g/kg
!                            TYPE:       REAL(fp)
!                            DIMENSION:  Scalar or any rank
!                            ATTRIBUTES: INTENT(IN)
!
!       Mixing_Ratio_TL:     Tangent-linear mass mixing ratio of gas.
!                            UNITS:      g/kg
!                            TYPE:       REAL(fp)
!                            DIMENSION:  Same as Mixing_Ratio argument.
!                            ATTRIBUTES: INTENT(IN)
!
! OUTPUTS:
!       Specific_Amount_TL:  Tangent-linear specific amount of gas.
!                            Set to 0.0 if input mixing ratio < 0.0
!                            UNITS:      g/kg
!                            TYPE:       REAL(fp)
!                            DIMENSION:  Same as Mixing_Ratio argument
!                            ATTRIBUTES: INTENT(OUT)
!
! OPTIONAL INPUT:
!       Water_Vapor:         Water vapor specific humidity. If this argument is
!                            not supplied, the input MIXING_RATIO argument is
!                            assumed to be water vapor.
!                            Output set to 0.0 if < 0.0
!                            UNITS:      g/kg
!                            TYPE:       REAL(fp)
!                            DIMENSION:  Same as Mixing_Ratio argument
!                            ATTRIBUTES: OPTIONAL, INTENT(IN)
!
!       Water_Vapor_TL:      Tangent-linear specific humidty. This argument is
!                            ignored if the optional Water_Vapor argument is
!                            not supplied also.
!                            UNITS:      g/kg
!                            TYPE:       REAL(fp)
!                            DIMENSION:  Same as Mixing_Ratio argument
!                            ATTRIBUTES: OPTIONAL, INTENT(IN)
!
! CREATION HISTORY:
!       Written by:     Paul van Delst, 14-Aug-2009
!                       paul.vandelst@noaa.gov
!:sdoc-:
!------------------------------------------------------------------------------

  ELEMENTAL SUBROUTINE MR_to_SA_TL( &
    Mixing_Ratio      , &  ! FWD Input
    Mixing_Ratio_TL   , &  ! TL Input
    Specific_Amount_TL, &  ! TL Output
    Water_Vapor       , &  ! Optional FWD Input
    Water_Vapor_TL      )  ! Optional TL  Input
    ! Arguments
    REAL(fp),           INTENT(IN)  :: Mixing_Ratio
    REAL(fp),           INTENT(IN)  :: Mixing_Ratio_TL
    REAL(fp),           INTENT(OUT) :: Specific_Amount_TL
    REAL(fp), OPTIONAL, INTENT(IN)  :: Water_Vapor
    REAL(fp), OPTIONAL, INTENT(IN)  :: Water_Vapor_TL

    ! Error checks
    ! ...Zero output for -ve input
    IF ( Mixing_Ratio < ZERO ) THEN
      Specific_Amount_TL = ZERO
      RETURN
    ENDIF

    ! Tangent-linear conversion
    IF ( PRESENT(Water_Vapor) ) THEN
      ! ...Zero output for -ve input
      IF ( Water_Vapor < ZERO ) THEN
        Specific_Amount_TL = ZERO
        RETURN
      ENDIF
      ! ...Gas other than water vapour
      Specific_Amount_TL = ((ONE - (G_TO_KG * Water_Vapor)) * Mixing_Ratio_TL)
      IF ( PRESENT(Water_Vapor_TL) ) THEN
        Specific_Amount_TL = Specific_Amount_TL - (G_TO_KG * Mixing_Ratio * Water_Vapor_TL)
      END IF
    ELSE
      ! ...Water vapour
      Specific_Amount_TL = Mixing_Ratio_TL / (ONE + (G_TO_KG * Mixing_Ratio))**2
    END IF

  END SUBROUTINE MR_to_SA_TL

!------------------------------------------------------------------------------
!:sdoc+:
! NAME:
!       MR_to_SA_AD
!
! PURPOSE:
!       Adjoint form of elemental subroutine to convert gas concentrations
!       from mass mixing ratio in g/kg to specific amount in g/kg.
!
! CALLING SEQUENCE:
!       CALL MR_to_SA_AD( Mixing_Ratio                   , &  ! FWD Input
!                         Specific_Amount_AD             , &  ! AD  Input
!                         Mixing_Ratio_AD                , &  ! AD  Output
!                         Water_Vapor    = Water_Vapor   , &  ! Optional FWD input
!                         Water_Vapor_AD = Water_Vapor_AD  )  ! Optional AD  output
!
! INPUTS:
!       Mixing_Ratio:        Mass mixing ratio of gas.
!                            Must be > or = 0.0
!                            UNITS:      g/kg
!                            TYPE:       REAL(fp)
!                            DIMENSION:  Scalar or any rank
!                            ATTRIBUTES: INTENT(IN)
!
!       Specific_Amount_AD:  Adjoint specific amount of gas.
!                            *** SET TO ZERO ON EXIT ***
!                            UNITS:      g/kg
!                            TYPE:       REAL(fp)
!                            DIMENSION:  Same as Mixing_Ratio argument
!                            ATTRIBUTES: INTENT(IN OUT)
!
! OUTPUTS:
!       Mixing_Ratio_AD:     Adjoint mass mixing ratio of gas.
!                            *** MUST HAVE VALUE ON ENTRY ***
!                            UNITS:      g/kg
!                            TYPE:       REAL(fp)
!                            DIMENSION:  Same as Mixing_Ratio argument.
!                            ATTRIBUTES: INTENT(IN OUT)
!
! OPTIONAL INPUT:
!       Water_Vapor:         Water vapor specific humidity. If this argument is
!                            not supplied, the input MIXING_RATIO argument is
!                            assumed to be water vapor.
!                            Output set to 0.0 if < 0.0
!                            UNITS:      g/kg
!                            TYPE:       REAL(fp)
!                            DIMENSION:  Same as Mixing_Ratio argument
!                            ATTRIBUTES: OPTIONAL, INTENT(IN)
!
! OPTIONAL OUTPUT:
!       Water_Vapor_AD:      Adjoint specific humidty. This argument is
!                            ignored if the optional Water_Vapor argument is
!                            not supplied also.
!                            *** MUST HAVE VALUE ON ENTRY ***
!                            UNITS:      g/kg
!                            TYPE:       REAL(fp)
!                            DIMENSION:  Same as Mixing_Ratio argument
!                            ATTRIBUTES: OPTIONAL, INTENT(IN OUT)
!
! CREATION HISTORY:
!       Written by:     Paul van Delst, 14-Aug-2009
!                       paul.vandelst@noaa.gov
!:sdoc-:
!------------------------------------------------------------------------------

  ELEMENTAL SUBROUTINE MR_to_SA_AD( &
    Mixing_Ratio      , &  ! FWD Input
    Specific_Amount_AD, &  ! AD Input
    Mixing_Ratio_AD   , &  ! AD Output
    Water_Vapor       , &  ! Optional FWD Input
    Water_Vapor_AD      )  ! Optional AD  Output
    ! Arguments
    REAL(fp),           INTENT(IN)     :: Mixing_Ratio
    REAL(fp),           INTENT(IN OUT) :: Specific_Amount_AD
    REAL(fp),           INTENT(IN OUT) :: Mixing_Ratio_AD
    REAL(fp), OPTIONAL, INTENT(IN)     :: Water_Vapor
    REAL(fp), OPTIONAL, INTENT(IN OUT) :: Water_Vapor_AD

    ! Error checks
    ! ...Zero output for -ve input
    IF ( Mixing_Ratio < ZERO ) THEN
      Specific_Amount_AD = ZERO
      Mixing_Ratio_AD    = ZERO
      RETURN
    ENDIF

    ! Tangent-linear conversion
    IF ( PRESENT(Water_Vapor) .AND. PRESENT(Water_Vapor_AD) ) THEN
      ! ...Zero output for -ve input
      IF ( Water_Vapor < ZERO ) THEN
        Specific_Amount_AD = ZERO
        Mixing_Ratio_AD    = ZERO
        Water_Vapor_AD     = ZERO
        RETURN
      ENDIF
      ! ...Gas other than water vapour
      Water_Vapor_AD  = Water_Vapor_AD - (G_TO_KG * Mixing_Ratio * Specific_Amount_AD)
      Mixing_Ratio_AD = Mixing_Ratio_AD + ((ONE - (G_TO_KG * Water_Vapor)) * Specific_Amount_AD)
    ELSE
      ! ...Water vapour
      Mixing_Ratio_AD = Mixing_Ratio_AD + (Specific_Amount_AD / (ONE + (G_TO_KG * Mixing_Ratio))**2)
    END IF
    Specific_Amount_AD = ZERO

  END SUBROUTINE MR_to_SA_AD

END MODULE SA_MR

max_stars_repo_name: hsbadr/crtm

max_issues_repo_name: hsbadr/crtm

max_forks_repo_name: hsbadr/crtm


! $Id$

module write_netcdf

contains

  ! Atmos test write routine for 'knowngood's
  subroutine buildtest_atm_write()

    use genie_control, only: outputdir_name,ilon1_atm,ilat1_atm,BUFSIZ
    use genie_global,  only: surf_tstarinst_atm,netlong_atm,alon1_atm, &
         alat1_atm,koverall

    ! locals
    character(len=BUFSIZ)                          :: fname
    character(len=BUFSIZ)                          :: koverallstring
    character(len=BUFSIZ),allocatable,dimension(:) :: varnames
    real,allocatable,dimension(:,:,:)              :: data

    ! construct filename
    write(koverallstring,'(i6.6)') koverall
    fname=trim(outputdir_name)//'/genie_buildtest_atm_' &
         //trim(koverallstring)//'_.nc'

    allocate(data(ilon1_atm,ilat1_atm,2))
    allocate(varnames(2))
    data(:,:,1) = surf_tstarinst_atm
    varnames(1) = 'surf_tstarinst_atm'
    data(:,:,2) = netlong_atm
    varnames(2) = 'netlong_atm'

    call write_netcdf_genie(fname,data,varnames,alon1_atm,alat1_atm)
    deallocate(data)
    deallocate(varnames)

  end subroutine buildtest_atm_write

  ! Ocean test write routine for 'knowngood's
  subroutine buildtest_ocn_write()

    use genie_control, only: outputdir_name,ilon1_ocn,ilat1_ocn,BUFSIZ
    use genie_global,  only: tstar_ocn,netlong_ocn,alon1_ocn,alat1_ocn, &
         koverall

    ! locals
    character(len=BUFSIZ)                          :: fname
    character(len=BUFSIZ)                          :: koverallstring
    character(len=BUFSIZ),allocatable,dimension(:) :: varnames
    real,allocatable,dimension(:,:,:)              :: data

    ! construct filename
    write(koverallstring,'(i6.6)') koverall
    fname=trim(outputdir_name)//'/genie_buildtest_ocn_' &
         //trim(koverallstring)//'_.nc'

    allocate(data(ilon1_ocn,ilat1_ocn,2))
    allocate(varnames(2))
    data(:,:,1) = tstar_ocn
    varnames(1) = 'tstar_ocn'
    data(:,:,2) = netlong_ocn
    varnames(2) = 'netlong_ocn'

    call write_netcdf_genie(fname,data,varnames,alon1_ocn,alat1_ocn)
    deallocate(data)
    deallocate(varnames)

  end subroutine buildtest_ocn_write

  subroutine write_netcdf_genie(filename,data1_in,invarnames,alon1_in,alat1_in)

!     ***********************************************************
!     djl, 14/11/2003
!     subroutine to enable instantaneous easy output of
!     grads/ferret-friendly single-level netcdf data
!
!     this version for the genie grid (ocean or atmosphere).
!     (see also write_netcdf_igcm.f and write_netcdf_std.f in genie-igcm3)
!
!     [filename    = output file name]
!     [data1       = data to be output]
!     [invarnames  = name to be given to the netcdf variable]
!     [alon1,alat1 = values of longitudes and latitudes]
!     [nx,ny       = number of longitudes and latitudes]
!
!     dagw, 17/11/2004
!     made into an f90 module
!     ***********************************************************

!!!    use precision
    use genie_control

    implicit none

    include 'netcdf.inc'

    ! locals
    integer                                           :: nx,ny  !nx-lon,ny-lat
    integer                                           :: ndim,nvar
    integer, dimension(nall)                          :: natts,nattsvar,vdims,ndims
    integer, dimension(nmaxdims,nall)                 :: vadims
    character(len=BUFSIZ), dimension(nall,nfiles)     :: dimname,varname
    character(len=BUFSIZ), dimension(2,nmaxdims,nall) :: attdimname,attvarname
    integer, dimension(nall,nfiles)                   :: iddim,idvar

    !     GENIE GRID:
    real, intent(in), dimension(:)             :: alon1_in
    real, intent(in), dimension(:)             :: alat1_in
    real, allocatable, dimension(:) :: alon1,alat1

    !     DATA:
    real, intent(in), dimension(:,:,:)             :: data1_in
    real, allocatable, dimension(:,:,:) :: data1

    character(len=*), intent(in)               :: filename
    character(len=*), intent(in), dimension(:) :: invarnames

    !     NETCDF + AWI STUFF:
    integer :: ncid,loc_dim

    !     LOOPING:
    integer :: i,j,k

    integer, dimension(3) :: dataShape
    integer, dimension(1) :: varNamesShape,lonShape,latShape

    ! determine extent of longititude and latitude arrays
    lonShape = shape(alon1_in)
    latShape = shape(alat1_in)
    nx = lonShape(1)
    ny = latShape(1)

    ! determine extent of data and varaible name arrays
    dataShape = shape(data1_in)
    varNamesShape = shape(invarnames)

    ! test for consistency
    ! NB **ASSUME** THAT DATA ARRAY IS (LON:LAT:NUMVARS)

    if (dataShape(1) .ne. nx) then
       print*,"Error:",__LINE__,__FILE__,"Longitude inconsistency"
       stop
    else if (dataShape(2) .ne. ny) then
       print*,"Error:",__LINE__,__FILE__,"Latitude inconsistency"
       stop
    else if (dataShape(3) .ne. varNamesShape(1)) then
       print*,"Error:",__LINE__,__FILE__,"Num vars inconsistency"
       stop
    end if

    nvar = varNamesShape(1)

    ! allocate internal storage
    allocate(alon1(nx))
    allocate(alat1(ny))
    allocate(data1(nx,ny,nvar))

    ! copy values--reasons of floating pt resolution
    do k=1,nvar
       do j=1,ny
          do i=1,nx
             data1(i,j,k)=data1_in(i,j,k)
          enddo
       enddo
    enddo

    do i=1,nx
       alon1(i)=alon1_in(i)
    enddo

    do j=1,ny
       alat1(j)=alat1_in(j)
    enddo

    ! Set up dimension information
    ! NB **ASSUME** DATA ON SAME LON-LAT GRID ONLY
    ! COPIED FROM INITIALISE_ATMOS.F
    ndim=2
    dimname(1,1)='longitude'
    ndims(1)=nx
    natts(1)=2
    attdimname(1,1,1)='long_name'
    attdimname(2,1,1)='longitude'
    attdimname(1,2,1)='units'
    attdimname(2,2,1)='degrees east'

    dimname(2,1)='latitude'
    ndims(2)=ny
    natts(2)=2
    attdimname(1,1,2)='long_name'
    attdimname(2,1,2)='latitude'
    attdimname(1,2,2)='units'
    attdimname(2,2,2)='degrees north'

    ! Setup variable information
    do k=1,nvar
       varname(k,1)=trim(invarnames(k))
       vdims(k)=2
       vadims(1,k)=loc_dim('longitude',dimname,nall)
       vadims(2,k)=loc_dim('latitude',dimname,nall)
       nattsvar(k)=2
       attvarname(1,1,k)='long_name'
       attvarname(2,1,k)='data from the model'
       attvarname(1,2,k)='units'
       attvarname(2,2,k)='no units'
    enddo

    call ininc(trim(filename),       &
         nmaxdims,ndim,nvar,         &
         natts,nattsvar,             &
         vdims,vadims,ndims,         &
         dimname(1,1),varname(1,1),  &
         attdimname,attvarname,      &
         ncid,iddim(1,1),idvar(1,1))

    call writedim(ncid,iddim(1,1),alon1)
    call writedim(ncid,iddim(2,1),alat1)

    do k=1,nvar
       call writevar(ncid,                                 &
            idvar(loc_dim(invarnames(k),varname(1,1),nall),1), &
            data1(:,:,k) )
    enddo

    call closenc(ncid)

    deallocate(alon1)
    deallocate(alat1)
    deallocate(data1)

    return
  end subroutine write_netcdf_genie
end module write_netcdf


text:
      SUBROUTINE HWSSSP (TS,TF,M,MBDCND,BDTS,BDTF,PS,PF,N,NBDCND,BDPS,
     1                   BDPF,ELMBDA,F,IDIMF,PERTRB,IERROR,W)
C
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     *                                                               *
C     *                        F I S H P A K                          *
C     *                                                               *
C     *                                                               *
C     *     A PACKAGE OF FORTRAN SUBPROGRAMS FOR THE SOLUTION OF      *
C     *                                                               *
C     *      SEPARABLE ELLIPTIC PARTIAL DIFFERENTIAL EQUATIONS        *
C     *                                                               *
C     *                  (VERSION 3.1 , OCTOBER 1980)                  *
C     *                                                               *
C     *                             BY                                *
C     *                                                               *
C     *        JOHN ADAMS, PAUL SWARZTRAUBER AND ROLAND SWEET         *
C     *                                                               *
C     *                             OF                                *
C     *                                                               *
C     *         THE NATIONAL CENTER FOR ATMOSPHERIC RESEARCH          *
C     *                                                               *
C     *                BOULDER, COLORADO  (80307)  U.S.A.             *
C     *                                                               *
C     *                   WHICH IS SPONSORED BY                       *
C     *                                                               *
C     *              THE NATIONAL SCIENCE FOUNDATION                  *
C     *                                                               *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C     SUBROUTINE HWSSSP SOLVES A FINITE DIFFERENCE APPROXIMATION TO THE
C     HELMHOLTZ EQUATION IN SPHERICAL COORDINATES AND ON THE SURFACE OF
C     THE UNIT SPHERE (RADIUS OF 1):
C
C          (1/SIN(THETA))(D/DTHETA)(SIN(THETA)(DU/DTHETA))
C
C             + (1/SIN(THETA)**2)(D/DPHI)(DU/DPHI)
C
C             + LAMBDA*U = F(THETA,PHI)
C
C     WHERE THETA IS COLATITUDE AND PHI IS LONGITUDE.
C
C     * * * * * * * *    PARAMETER DESCRIPTION     * * * * * * * * * *
C
C             * * * * * *   ON INPUT    * * * * * *
C
C     TS,TF
C       THE RANGE OF THETA (COLATITUDE), I.E., TS .LE. THETA .LE. TF.
C       TS MUST BE LESS THAN TF.  TS AND TF ARE IN RADIANS.  A TS OF
C       ZERO CORRESPONDS TO THE NORTH POLE AND A TF OF PI CORRESPONDS TO
C       THE SOUTH POLE.
C
C     * * * * * * * * * * * * * * IMPORTANT * * * * * * * * * * * * * *
C
C     IF TF IS EQUAL TO PI THEN IT MUST BE COMPUTED USING THE STATEMENT
C     TF = PIMACH(DUM). THIS INSURES THAT TF IN THE USERS PROGRAM IS
C     EQUAL TO PI IN THIS PROGRAM WHICH PERMITS SEVERAL TESTS OF THE
C     INPUT PARAMETERS THAT OTHERWISE WOULD NOT BE POSSIBLE.
C
C
C     M
C       THE NUMBER OF PANELS INTO WHICH THE INTERVAL (TS,TF) IS
C       SUBDIVIDED.  HENCE, THERE WILL BE M+1 GRID POINTS IN THE
C       THETA-DIRECTION GIVEN BY THETA(I) = (I-1)DTHETA+TS FOR
C       I = 1,2,...,M+1, WHERE DTHETA = (TF-TS)/M IS THE PANEL WIDTH.
C       M MUST BE GREATER THAN 5
C
C     MBDCND
C       INDICATES THE TYPE OF BOUNDARY CONDITION AT THETA = TS AND
C       THETA = TF.
C
C       = 1  IF THE SOLUTION IS SPECIFIED AT THETA = TS AND THETA = TF.
C       = 2  IF THE SOLUTION IS SPECIFIED AT THETA = TS AND THE
C            DERIVATIVE OF THE SOLUTION WITH RESPECT TO THETA IS
C            SPECIFIED AT THETA = TF (SEE NOTE 2 BELOW).
C       = 3  IF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO THETA IS
C            SPECIFIED AT THETA = TS AND THETA = TF (SEE NOTES 1,2
C            BELOW).
C       = 4  IF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO THETA IS
C            SPECIFIED AT THETA = TS (SEE NOTE 1 BELOW) AND THE
C            SOLUTION IS SPECIFIED AT THETA = TF.
C       = 5  IF THE SOLUTION IS UNSPECIFIED AT THETA = TS = 0 AND THE
C            SOLUTION IS SPECIFIED AT THETA = TF.
C       = 6  IF THE SOLUTION IS UNSPECIFIED AT THETA = TS = 0 AND THE
C            DERIVATIVE OF THE SOLUTION WITH RESPECT TO THETA IS
C            SPECIFIED AT THETA = TF (SEE NOTE 2 BELOW).
C       = 7  IF THE SOLUTION IS SPECIFIED AT THETA = TS AND THE
C            SOLUTION IS UNSPECIFIED AT THETA = TF = PI.
C       = 8  IF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO THETA IS
C            SPECIFIED AT THETA = TS (SEE NOTE 1 BELOW) AND THE
C            SOLUTION IS UNSPECIFIED AT THETA = TF = PI.
C       = 9  IF THE SOLUTION IS UNSPECIFIED AT THETA = TS = 0 AND
C            THETA = TF = PI.
C
C       NOTES:  1.  IF TS = 0, DO NOT USE MBDCND = 3,4, OR 8, BUT
C                   INSTEAD USE MBDCND = 5,6, OR 9  .
C               2.  IF TF = PI, DO NOT USE MBDCND = 2,3, OR 6, BUT
C                   INSTEAD USE MBDCND = 7,8, OR 9  .
C
C     BDTS
C       A ONE-DIMENSIONAL ARRAY OF LENGTH N+1 THAT SPECIFIES THE VALUES
C       OF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO THETA AT
C       THETA = TS.  WHEN MBDCND = 3,4, OR 8,
C
C            BDTS(J) = (D/DTHETA)U(TS,PHI(J)), J = 1,2,...,N+1  .
C
C       WHEN MBDCND HAS ANY OTHER VALUE, BDTS IS A DUMMY VARIABLE.
C
C     BDTF
C       A ONE-DIMENSIONAL ARRAY OF LENGTH N+1 THAT SPECIFIES THE VALUES
C       OF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO THETA AT
C       THETA = TF.  WHEN MBDCND = 2,3, OR 6,
C
C            BDTF(J) = (D/DTHETA)U(TF,PHI(J)), J = 1,2,...,N+1  .
C
C       WHEN MBDCND HAS ANY OTHER VALUE, BDTF IS A DUMMY VARIABLE.
C
C     PS,PF
C       THE RANGE OF PHI (LONGITUDE), I.E., PS .LE. PHI .LE. PF.  PS
C       MUST BE LESS THAN PF.  PS AND PF ARE IN RADIANS.  IF PS = 0 AND
C       PF = 2*PI, PERIODIC BOUNDARY CONDITIONS ARE USUALLY PRESCRIBED.
C
C     * * * * * * * * * * * * * * IMPORTANT * * * * * * * * * * * * * *
C
C     IF PF IS EQUAL TO 2*PI THEN IT MUST BE COMPUTED USING THE
C     STATEMENT PF = 2.*PIMACH(DUM). THIS INSURES THAT PF IN THE USERS
C     PROGRAM IS EQUAL TO 2*PI IN THIS PROGRAM WHICH PERMITS TESTS OF
C     THE INPUT PARAMETERS THAT OTHERWISE WOULD NOT BE POSSIBLE.
C
C
C     N
C       THE NUMBER OF PANELS INTO WHICH THE INTERVAL (PS,PF) IS
C       SUBDIVIDED.  HENCE, THERE WILL BE N+1 GRID POINTS IN THE
C       PHI-DIRECTION GIVEN BY PHI(J) = (J-1)DPHI+PS  FOR
C       J = 1,2,...,N+1, WHERE DPHI = (PF-PS)/N IS THE PANEL WIDTH.
C       N MUST BE GREATER THAN 4
C
C     NBDCND
C       INDICATES THE TYPE OF BOUNDARY CONDITION AT PHI = PS AND
C       PHI = PF.
C
C       = 0  IF THE SOLUTION IS PERIODIC IN PHI, I.E.,
C            U(I,J) = U(I,N+J).
C       = 1  IF THE SOLUTION IS SPECIFIED AT PHI = PS AND PHI = PF
C            (SEE NOTE BELOW).
C       = 2  IF THE SOLUTION IS SPECIFIED AT PHI = PS (SEE NOTE BELOW)
C            AND THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO PHI IS
C            SPECIFIED AT PHI = PF.
C       = 3  IF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO PHI IS
C            SPECIFIED AT PHI = PS AND PHI = PF.
C       = 4  IF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO PHI IS
C            SPECIFIED AT PS AND THE SOLUTION IS SPECIFIED AT PHI = PF
C            (SEE NOTE BELOW).
C
C       NOTE:  NBDCND = 1,2, OR 4 CANNOT BE USED WITH
C              MBDCND = 5,6,7,8, OR 9 (THE FORMER INDICATES THAT THE
C                       SOLUTION IS SPECIFIED AT A POLE, THE LATTER
C                       INDICATES THAT THE SOLUTION IS UNSPECIFIED).
C                       USE INSTEAD
C              MBDCND = 1 OR 2  .
C
C     BDPS
C       A ONE-DIMENSIONAL ARRAY OF LENGTH M+1 THAT SPECIFIES THE VALUES
C       OF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO PHI AT
C       PHI = PS.  WHEN NBDCND = 3 OR 4,
C
C            BDPS(I) = (D/DPHI)U(THETA(I),PS), I = 1,2,...,M+1  .
C
C       WHEN NBDCND HAS ANY OTHER VALUE, BDPS IS A DUMMY VARIABLE.
C
C     BDPF
C       A ONE-DIMENSIONAL ARRAY OF LENGTH M+1 THAT SPECIFIES THE VALUES
C       OF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO PHI AT
C       PHI = PF.  WHEN NBDCND = 2 OR 3,
C
C            BDPF(I) = (D/DPHI)U(THETA(I),PF), I = 1,2,...,M+1  .
C
C       WHEN NBDCND HAS ANY OTHER VALUE, BDPF IS A DUMMY VARIABLE.
C
C     ELMBDA
C       THE CONSTANT LAMBDA IN THE HELMHOLTZ EQUATION.  IF
C       LAMBDA .GT. 0, A SOLUTION MAY NOT EXIST.  HOWEVER, HWSSSP WILL
C       ATTEMPT TO FIND A SOLUTION.
C
C     F
C       A TWO-DIMENSIONAL ARRAY THAT SPECIFIES THE VALUE OF THE RIGHT
C       SIDE OF THE HELMHOLTZ EQUATION AND BOUNDARY VALUES (IF ANY).
C       FOR I = 2,3,...,M  AND  J = 2,3,...,N
C
C            F(I,J) = F(THETA(I),PHI(J)).
C
C       ON THE BOUNDARIES F IS DEFINED BY
C
C            MBDCND   F(1,J)            F(M+1,J)
C            ------   ------------      ------------
C
C              1      U(TS,PHI(J))      U(TF,PHI(J))
C              2      U(TS,PHI(J))      F(TF,PHI(J))
C              3      F(TS,PHI(J))      F(TF,PHI(J))
C              4      F(TS,PHI(J))      U(TF,PHI(J))
C              5      F(0,PS)           U(TF,PHI(J))   J = 1,2,...,N+1
C              6      F(0,PS)           F(TF,PHI(J))
C              7      U(TS,PHI(J))      F(PI,PS)
C              8      F(TS,PHI(J))      F(PI,PS)
C              9      F(0,PS)           F(PI,PS)
C
C            NBDCND   F(I,1)            F(I,N+1)
C            ------   --------------    --------------
C
C              0      F(THETA(I),PS)    F(THETA(I),PS)
C              1      U(THETA(I),PS)    U(THETA(I),PF)
C              2      U(THETA(I),PS)    F(THETA(I),PF)   I = 1,2,...,M+1
C              3      F(THETA(I),PS)    F(THETA(I),PF)
C              4      F(THETA(I),PS)    U(THETA(I),PF)
C
C       F MUST BE DIMENSIONED AT LEAST (M+1)*(N+1).
C       NOTE
C
C       IF THE TABLE CALLS FOR BOTH THE SOLUTION U AND THE RIGHT SIDE F
C       AT  A CORNER THEN THE SOLUTION MUST BE SPECIFIED.
C
C
C     IDIMF
C       THE ROW (OR FIRST) DIMENSION OF THE ARRAY F AS IT APPEARS IN THE
C       PROGRAM CALLING HWSSSP.  THIS PARAMETER IS USED TO SPECIFY THE
C       VARIABLE DIMENSION OF F.  IDIMF MUST BE AT LEAST M+1  .
C
C     W
C       A ONE-DIMENSIONAL ARRAY THAT MUST BE PROVIDED BY THE USER FOR
C       WORK SPACE. W MAY REQUIRE UP TO 4*(N+1)+(16+INT(LOG2(N+1)))(M+1)
C       LOCATIONS. THE ACTUAL NUMBER OF LOCATIONS USED IS COMPUTED BY
C       HWSSSP AND IS OUTPUT IN LOCATION W(1). INT( ) DENOTES THE
C       FORTRAN INTEGER FUNCTION.
C
C
C     * * * * * * * * * *     ON OUTPUT     * * * * * * * * * *
C
C     F
C       CONTAINS THE SOLUTION U(I,J) OF THE FINITE DIFFERENCE
C       APPROXIMATION FOR THE GRID POINT (THETA(I),PHI(J)),
C       I = 1,2,...,M+1,   J = 1,2,...,N+1  .
C
C     PERTRB
C       IF ONE SPECIFIES A COMBINATION OF PERIODIC, DERIVATIVE OR
C       UNSPECIFIED BOUNDARY CONDITIONS FOR A POISSON EQUATION
C       (LAMBDA = 0), A SOLUTION MAY NOT EXIST.  PERTRB IS A CONSTANT,
C       CALCULATED AND SUBTRACTED FROM F, WHICH ENSURES THAT A SOLUTION
C       EXISTS.  HWSSSP THEN COMPUTES THIS SOLUTION, WHICH IS A LEAST
C       SQUARES SOLUTION TO THE ORIGINAL APPROXIMATION.  THIS SOLUTION
C       IS NOT UNIQUE AND IS UNNORMALIZED. THE VALUE OF PERTRB SHOULD
C       BE SMALL COMPARED TO THE RIGHT SIDE F. OTHERWISE , A SOLUTION
C       IS OBTAINED TO AN ESSENTIALLY DIFFERENT PROBLEM. THIS COMPARISON
C       SHOULD ALWAYS BE MADE TO INSURE THAT A MEANINGFUL SOLUTION HAS
C       BEEN OBTAINED
C
C     IERROR
C       AN ERROR FLAG THAT INDICATES INVALID INPUT PARAMETERS.  EXCEPT
C       FOR NUMBERS 0 AND 8, A SOLUTION IS NOT ATTEMPTED.
C
C       = 0  NO ERROR
C       = 1  TS.LT.0 OR TF.GT.PI
C       = 2  TS.GE.TF
C       = 3  MBDCND.LT.1 OR MBDCND.GT.9
C       = 4  PS.LT.0 OR PS.GT.PI+PI
C       = 5  PS.GE.PF
C       = 6  N.LT.5
C       = 7  M.LT.5
C       = 8  NBDCND.LT.0 OR NBDCND.GT.4
C       = 9  ELMBDA.GT.0
C       = 10 IDIMF.LT.M+1
C       = 11 NBDCND EQUALS 1,2 OR 4 AND MBDCND.GE.5
C       = 12 TS.EQ.0 AND MBDCND EQUALS 3,4 OR 8
C       = 13 TF.EQ.PI AND MBDCND EQUALS 2,3 OR 6
C       = 14 MBDCND EQUALS 5,6 OR 9 AND TS.NE.0
C       = 15 MBDCND.GE.7 AND TF.NE.PI
C
C       SINCE THIS IS THE ONLY MEANS OF INDICATING A POSSIBLY INCORRECT
C       CALL TO HWSSSP, THE USER SHOULD TEST IERROR AFTER A CALL.
C
C     W
C       CONTAINS INTERMEDIATE VALUES THAT MUST NOT BE DESTROYED IF
C       HWSSSP WILL BE CALLED AGAIN WITH INTL = 1. W(1) CONTAINS THE
C       REQUIRED LENGTH OF W .
C
C
C
C     * * * * * * *   PROGRAM SPECIFICATIONS    * * * * * * * * * * * *
C
C     DIMENSION OF   BDTS(N+1),BDTF(N+1),BDPS(M+1),BDPF(M+1),
C     ARGUMENTS      F(IDIMF,N+1),W(SEE ARGUMENT LIST)
C
C     LATEST         JANUARY 1978
C     REVISION
C
C
C     SUBPROGRAMS    HWSSSP,HWSSS1,GENBUN,POISD2,POISN2,POISP2,COSGEN,ME
C     REQUIRED       TRIX,TRI3,PIMACH
C
C     SPECIAL        NONE
C     CONDITIONS
C
C     COMMON         NONE
C     BLOCKS
C
C     I/O            NONE
C
C     PRECISION      SINGLE
C
C     SPECIALIST     PAUL SWARZTRAUBER
C
C     LANGUAGE       FORTRAN
C
C     HISTORY        VERSION 1 - SEPTEMBER 1973
C                    VERSION 2 - APRIL     1976
C                    VERSION 3 - JANUARY   1978
C
C     ALGORITHM      THE ROUTINE DEFINES THE FINITE DIFFERENCE
C                    EQUATIONS, INCORPORATES BOUNDARY DATA, AND ADJUSTS
C                    THE RIGHT SIDE OF SINGULAR SYSTEMS AND THEN CALLS
C                    GENBUN TO SOLVE THE SYSTEM.
C
C     SPACE
C     REQUIRED       CONTROL DATA 7600
C
C     TIMING AND        THE EXECUTION TIME T ON THE NCAR CONTROL DATA
C     ACCURACY       7600 FOR SUBROUTINE HWSSSP IS ROUGHLY PROPORTIONAL
C                    TO M*N*LOG2(N), BUT ALSO DEPENDS ON THE INPUT
C                    PARAMETERS NBDCND AND MBDCND.  SOME TYPICAL VALUES
C                    ARE LISTED IN THE TABLE BELOW.
C                       THE SOLUTION PROCESS EMPLOYED RESULTS IN A LOSS
C                    OF NO MORE THAN THREE SIGNIFICANT DIGITS FOR N AND
C                    M AS LARGE AS 64.  MORE DETAILED INFORMATION ABOUT
C                    ACCURACY CAN BE FOUND IN THE DOCUMENTATION FOR
C                    SUBROUTINE GENBUN WHICH IS THE ROUTINE THAT
C                    SOLVES THE FINITE DIFFERENCE EQUATIONS.
C
C
C                       M(=N)    MBDCND    NBDCND    T(MSECS)
C                       -----    ------    ------    --------
C
C                        32        0         0          31
C                        32        1         1          23
C                        32        3         3          36
C                        64        0         0         128
C                        64        1         1          96
C                        64        3         3         142
C
C     PORTABILITY    AMERICAN NATIONAL STANDARDS INSTITUTE FORTRAN.
C                    ALL MACHINE DEPENDENT CONSTANTS ARE LOCATED IN THE
C                    FUNCTION PIMACH.
C
C     REQUIRED       SIN,COS
C     RESIDENT
C     ROUTINES
C
C     REPERENCES     P. N. SWARZTRAUBER,THE DIRECT SOLUTION OF THE
C                    DISCRETE POISSON EQUATION ON THE SURFACE OF A
C                    SPHERE,S.I.A.M. J. NUMER. ANAL.,15(1974),PP 212-215
C
C                    SWARZTRAUBER,P. AND R. SWEET, 'EFFICIENT FORTRAN
C                    SUBPROGRAMS FOR THE SOLUTION OF ELLIPTIC EQUATIONS'
C                    NCAR TN/IA-109, JULY, 1975, 138 PP.
C
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C
      DIMENSION       F(IDIMF,1) ,BDTS(1)    ,BDTF(1)    ,BDPS(1)    ,
     1                BDPF(1)    ,W(1)
      NBR = NBDCND+1
      PI = PIMACH(DUM)
      TPI = 2.*PI
      IERROR = 0
      IF (TS.LT.0. .OR. TF.GT.PI) IERROR = 1
      IF (TS .GE. TF) IERROR = 2
      IF (MBDCND.LT.1 .OR. MBDCND.GT.9) IERROR = 3
      IF (PS.LT.0. .OR. PF.GT.TPI) IERROR = 4
      IF (PS .GE. PF) IERROR = 5
      IF (N .LT. 5) IERROR = 6
      IF (M .LT. 5) IERROR = 7
      IF (NBDCND.LT.0 .OR. NBDCND.GT.4) IERROR = 8
      IF (ELMBDA .GT. 0.) IERROR = 9
      IF (IDIMF .LT. M+1) IERROR = 10
      IF ((NBDCND.EQ.1 .OR. NBDCND.EQ.2 .OR. NBDCND.EQ.4) .AND.
     1    MBDCND.GE.5) IERROR = 11
      IF (TS.EQ.0. .AND.
     1    (MBDCND.EQ.3 .OR. MBDCND.EQ.4 .OR. MBDCND.EQ.8)) IERROR = 12
      IF (TF.EQ.PI .AND.
     1    (MBDCND.EQ.2 .OR. MBDCND.EQ.3 .OR. MBDCND.EQ.6)) IERROR = 13
      IF ((MBDCND.EQ.5 .OR. MBDCND.EQ.6 .OR. MBDCND.EQ.9) .AND.
     1    TS.NE.0.) IERROR = 14
      IF (MBDCND.GE.7 .AND. TF.NE.PI) IERROR = 15
      IF (IERROR.NE.0 .AND. IERROR.NE.9) RETURN
      CALL HWSSS1 (TS,TF,M,MBDCND,BDTS,BDTF,PS,PF,N,NBDCND,BDPS,BDPF,
     1             ELMBDA,F,IDIMF,PERTRB,W,W(M+2),W(2*M+3),W(3*M+4),
     2             W(4*M+5),W(5*M+6),W(6*M+7))
      W(1) = W(6*M+7)+FLOAT(6*(M+1))
      RETURN
      END

max_stars_repo_name: sbwhitecap/clocc-hg

max_issues_repo_name: sbwhitecap/clocc-hg

max_forks_repo_name: sbwhitecap/clocc-hg


text:
Freeborn Hall, formerly known as Memorial Union Assembly Hall, was UC Davis primary large event venue until the construction of the Mondavi Center. It was named for Stanley Freeborn Stanley B. Freeborn. Construction of the 1,800 seat building was completed in 1961. There are often conferences as well as the occasional speaker or concert. Nearby food: Coffee House. address(Located in, 38.542594, 121.750173) the basement, Lower Freeborn Hall, are such fine institutions such as The California Aggie, KDVS, Project Compost and AS Papers.

Due to seismic issues, Freeborn will be closed indefinitely in the Summer of 2014, and will be assessed by the Department of Campus Recreation and Unions to determine the scope and feasibility of renovation.

There is a tunnel connecting the MUGA to Lower Freeborn.

Freeborn has seen such prominent performers as Nine Inch Nails, The Grateful Dead, and former President Bill Clinton made a brief visit after his appearance at the Mondavi Center.

The basement was formerly home to University Haircutters.

 Photos

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search