Class radiation_two_stream_app
In: radiation/radiation_two_stream_app.f90

Methods

IDSchemeLongWave   IDSchemeShortWave   RadiationTwoStreamApp   RadiationTwoStreamApp   RadiationTwoStreamAppCore   RadiationTwoStreamAppWrapper   module_name   tridiag   tridiag1   tridiaginv   version  

Included Modules

dc_types gridset dc_message constants

Public Instance methods

RadiationTwoStreamApp( QeRatio, xyz_SSA, xyz_AF, SolarFluxTOA, xy_SurfAlbedo, IDScheme, xy_InAngle, xyr_OptDepBase, xyr_RadFlux )
Subroutine :
QeRatio :real(DP), intent(in )
xyz_SSA( 0:imax-1, 1:jmax, 1:kmax ) :real(DP), intent(in )
xyz_AF( 0:imax-1, 1:jmax, 1:kmax ) :real(DP), intent(in )
SolarFluxTOA :real(DP), intent(in )
xy_SurfAlbedo( 0:imax-1, 1:jmax ) :real(DP), intent(in )
IDScheme :integer , intent(in )
xy_InAngle(0:imax-1, 1:jmax) :real(DP), intent(in )
: sec (入射角). sec (angle of incidence)
xyr_OptDepBase( 0:imax-1, 1:jmax, 0:kmax ) :real(DP), intent(in )
xyr_RadFlux( 0:imax-1, 1:jmax, 0:kmax ) :real(DP), intent(out)

Alias for RadiationTwoStreamAppWrapper

RadiationTwoStreamApp( QeRatio, xyz_SSA, xyz_AF, SolarFluxTOA, xy_SurfAlbedo, IDScheme, xy_InAngle, xyr_OptDepBase, xyr_RadFlux, js, je )
Subroutine :
QeRatio :real(DP), intent(in )
xyz_SSA( 0:imax-1, 1:jmax, 1:kmax ) :real(DP), intent(in )
xyz_AF( 0:imax-1, 1:jmax, 1:kmax ) :real(DP), intent(in )
SolarFluxTOA :real(DP), intent(in )
xy_SurfAlbedo( 0:imax-1, 1:jmax ) :real(DP), intent(in )
IDScheme :integer , intent(in )
xy_InAngle(0:imax-1, 1:jmax) :real(DP), intent(in )
: sec (入射角). sec (angle of incidence)
xyr_OptDepBase( 0:imax-1, 1:jmax, 0:kmax ) :real(DP), intent(in )
xyr_RadFlux( 0:imax-1, 1:jmax, 0:kmax ) :real(DP), intent(out)
js :integer , intent(in )
je :integer , intent(in )

Alias for RadiationTwoStreamAppCore

Private Instance methods

IDSchemeLongWave
Constant :
IDSchemeLongWave = 2 :integer, parameter
IDSchemeShortWave
Constant :
IDSchemeShortWave = 1 :integer, parameter
RadiationTwoStreamAppCore( QeRatio, xyz_SSA, xyz_AF, SolarFluxTOA, xy_SurfAlbedo, IDScheme, xy_InAngle, xyr_OptDepBase, xyr_RadFlux, js, je )
Subroutine :
QeRatio :real(DP), intent(in )
xyz_SSA( 0:imax-1, 1:jmax, 1:kmax ) :real(DP), intent(in )
xyz_AF( 0:imax-1, 1:jmax, 1:kmax ) :real(DP), intent(in )
SolarFluxTOA :real(DP), intent(in )
xy_SurfAlbedo( 0:imax-1, 1:jmax ) :real(DP), intent(in )
IDScheme :integer , intent(in )
xy_InAngle(0:imax-1, 1:jmax) :real(DP), intent(in )
: sec (入射角). sec (angle of incidence)
xyr_OptDepBase( 0:imax-1, 1:jmax, 0:kmax ) :real(DP), intent(in )
xyr_RadFlux( 0:imax-1, 1:jmax, 0:kmax ) :real(DP), intent(out)
js :integer , intent(in )
je :integer , intent(in )

[Source]

  subroutine RadiationTwoStreamAppCore( QeRatio, xyz_SSA, xyz_AF, SolarFluxTOA, xy_SurfAlbedo, IDScheme, xy_InAngle, xyr_OptDepBase, xyr_RadFlux, js, je )

    ! USE statements
    !

    ! メッセージ出力
    ! Message output
    !
    use dc_message, only: MessageNotify

    !
    ! Physical constants settings
    !
    use constants, only: Grav, CpDry, PI, StB      ! $ \sigma_{SB} $ .
                                  ! ステファンボルツマン定数.
                                  ! Stefan-Boltzmann constant

    ! 格子点設定
    ! Grid points settings
    !
    use gridset, only: imax, jmax, kmax    ! 鉛直層数.
                               ! Number of vertical level

!!$  use pf_module , only : pfint_gq_array


    real(DP), intent(in ) :: QeRatio
    real(DP), intent(in ) :: xyz_SSA       ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP), intent(in ) :: xyz_AF        ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP), intent(in ) :: SolarFluxTOA
    real(DP), intent(in ) :: xy_SurfAlbedo ( 0:imax-1, 1:jmax )
    integer , intent(in ) :: IDScheme
    real(DP), intent(in ) :: xy_InAngle    (0:imax-1, 1:jmax)
                              ! sec (入射角).
                              ! sec (angle of incidence)
    real(DP), intent(in ) :: xyr_OptDepBase( 0:imax-1, 1:jmax, 0:kmax )
    real(DP), intent(out) :: xyr_RadFlux   ( 0:imax-1, 1:jmax, 0:kmax )

    integer , intent(in ) :: js
    integer , intent(in ) :: je


!!$    real(DP), intent(in ) :: gt         ( 0:imax-1, 1:jmax, 1:kmax )
!!$    real(DP), intent(in ) :: gts        ( 0:imax-1, 1:jmax )
!!$    real(DP), intent(in ) :: gph        ( 0:imax-1, 1:jmax, 0:kmax )

!!$    real(DP), intent(in ) :: emis  ( 0:imax-1, 1:jmax )
!!$    real(DP), intent(in ) :: wn1, wn2
!!$    integer , intent(in ) :: divnum


!!$  real(DP)    , intent(out) :: &
!!$    gor ( im, jm ), goru ( im, jm ), gord ( im, jm ), &
!!$    gsr ( im, jm ), gsru ( im, jm ), gsrd ( im, jm )

    !
    ! opdep       : Optical depth
    !
    real(DP) :: xyr_OptDep( 0:imax-1, 1:jmax, 0:kmax )

    !
    ! ssa2    : Delta-Function Adjusted Single-Scattering Albedo
    ! af2     : Delta-Function Adjusted Asymmetry Factor
    ! opdep2  : Delta-Function Adjusted Optical Depth
    !
    real(DP) :: xyz_SSA2  ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: xyz_AF2   ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: xyr_OpDep2( 0:imax-1, 1:jmax, 0:kmax )
    real(DP) :: xyr_Trans2( 0:imax-1, 1:jmax, 0:kmax )

    !
    ! gam?    : Coefficients of Generalized Radiative Transfer Equation
    !
    real(DP) :: xyz_Gam1( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: xyz_Gam2( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: xyz_Gam3( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: xyz_Gam4( 0:imax-1, 1:jmax, 1:kmax )

    !
    ! upfl      : Upward Flux
    ! dofl      : Downward Flux
    !
    real(DP) :: xyr_UwFlux( 0:imax-1, 1:jmax, 0:kmax )
    real(DP) :: xyr_DwFlux( 0:imax-1, 1:jmax, 0:kmax )

    !
    ! cosz     : cosine of solar zenith angle
    ! cosz2    : cosz squared
    !
    real(DP) :: xy_cosSZA     ( 0:imax-1, 1:jmax )
    real(DP) :: xy_cosSZAInv  ( 0:imax-1, 1:jmax )
    real(DP) :: xy_cosSZAInvsq( 0:imax-1, 1:jmax )

    !
    ! temporary variables
    !
    real(DP) :: xyz_DTau( 0:imax-1, 1:jmax, 1:kmax )

    integer  :: i, j, k, l
    integer  :: ms, me


    real(DP) :: xyz_Lambda ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: xyz_LSigma ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: xyaz_smalle( 0:imax-1, 1:jmax, 4, 1:kmax )

    !
    ! cupb      : upward C at bottom of layer
    ! cupt      : upward C at top of layer
    ! cdob      : downward C at bottom of layer
    ! cdot      : downward C at top of layer
    !
    real(DP) :: xyz_cupb( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: xyz_cupt( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: xyz_cdob( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: xyz_cdot( 0:imax-1, 1:jmax, 1:kmax )

    real(DP) :: aa_TridiagMtx1( 1:imax*jmax, 1:kmax*2 )
    real(DP) :: aa_TridiagMtx2( 1:imax*jmax, 1:kmax*2 )
    real(DP) :: aa_TridiagMtx3( 1:imax*jmax, 1:kmax*2 )
    real(DP) :: aa_Vec        ( 1:imax*jmax, 1:kmax*2 )

    real(DP) :: xy_tmpval( 0:imax-1, 1:jmax )

!!$  real(DP) :: mu1
!!$  real(DP) :: gth( im, jm, km+1 ), pfinth( im, jm, km+1 )
!!$  real(DP) :: b0( ijs:ije, 1, km ), b1( ijs:ije, 1, km )
!!$  real(DP) :: gemis


!!$    real(DP) :: inteup( 0:imax-1, 1:jmax, 0:kmax )
!!$    real(DP) :: intedo( 0:imax-1, 1:jmax, 0:kmax )
!!$    real(DP) :: tmpg, tmph, tmpj, tmpk, alp1, alp2, sig1, sig2
!!$    real(DP) :: tmpcos
!!$!      integer(i4b), parameter :: quadn = 3
!!$    integer(i4b), parameter :: quadn = 5
!!$!      integer(i4b), parameter :: quadn = 8
!!$    real(DP) :: qucos ( quadn ), qudcos( quadn )



!!$      call gauleg( 0.0d0, 1.0d0, quadn, qucos, qudcos )

    do k = 1, kmax
      do j = js, je
        do i = 0, imax-1
          if ( xyz_ssa(i,j,k) >= 1.0d0 ) then
            call MessageNotify( 'E', module_name, 'Single Scattering Albedo = %f', d = (/ xyz_ssa(i,j,k) /) )
          end if
        end do
      end do
    end do


    if( IDScheme .eq. IDSchemeShortWave ) then
      do j = js, je
        do i = 0, imax-1
          if ( xy_InAngle(i,j) > 0.0d0 ) then
            xy_cosSZA     (i,j) = 1.0d0 / xy_InAngle(i,j)
            xy_cosSZAInv  (i,j) = xy_InAngle(i,j)
            xy_cosSZAInvsq(i,j) = xy_cosSZAInv(i,j)**2
          else
            xy_cosSZA     (i,j) = 0.0d0
            xy_cosSZAInv  (i,j) = 0.0d0
            xy_cosSZAInvsq(i,j) = 0.0d0
          end if
        end do
      end do


!!$      gth(:,:,:) = 1.0d100

!!$    else
!!$
!!$      stop 'sw != 1 is not supported.'
!!$
!!$         do ij = ijs, ije
!!$            cosz ( ij, 1 ) = 1.0d100
!!$            cosz2( ij, 1 ) = 1.0d100
!!$         end do
!!$
!!$         do ij = ijs, ije
!!$            gth( ij, 1, 1 ) = gt( ij, 1, 1 )
!!$         end do
!!$         do k = 2, km+1-1
!!$            do ij = ijs, ije
!!$               gth( ij, 1, k ) = ( gt( ij, 1, k-1 ) + gt( ij, 1, k ) ) * 0.5d0
!!$            end do
!!$         end do
!!$         do ij = ijs, ije
!!$            gth( ij, 1, km+1 ) = gts( ij, 1 )
!!$         end do
!!$!         call pfint_array( wn1, wn2, divnum, im, jm, km+1, gth, pfinth, &
!!$!              ijs, ije )
!!$         call pfint_gq_array( wn1, wn2, divnum, im, jm, km+1, gth, pfinth, &
!!$              ijs, ije )

    else
      stop 'Unexpected sw number in twostreamapp.f90'
    end if


    !
    ! calculation of Dust Optical Depth
    !
    do k = 0, kmax
      do j = js, je
        xyr_OptDep(:,j,k) = xyr_OptDepBase(:,j,k) * QeRatio
      end do
    end do


    if( IDScheme .eq. IDSchemeShortWave ) then

      !
      ! Delta-Function Adjustment
      !
      do k = 1, kmax
        do j = js, je
          xyz_AF2 (:,j,k) = xyz_AF(:,j,k) / ( 1.0d0 + xyz_AF(:,j,k) )
          xyz_SSA2(:,j,k) =   ( 1.0d0 - xyz_AF(:,j,k)**2 ) * xyz_SSA(:,j,k) / ( 1.0d0 - xyz_SSA(:,j,k) * xyz_AF(:,j,k)**2 )
        end do
      end do

!!$      opdep2(:,:,:) = ( 1.0d0 - xyz_ssa * xyz_af**2 ) * opdep(:,:,:)
      do k = 0, kmax
        do j = js, je
          xyr_OpDep2(:,j,k) = 0.0d0
        end do
      end do
      do k = kmax-1, 0, -1
        do j = js, je
          xyr_OpDep2(:,j,k) = xyr_OpDep2(:,j,k+1) + ( 1.0d0 - xyz_SSA(:,j,k+1) * xyz_AF(:,j,k+1)**2 ) * ( xyr_OptDep(:,j,k) - xyr_OptDep(:,j,k+1) )
        end do
      end do

      do k = 0, kmax
        do j = js, je
          xyr_Trans2(:,j,k) = exp( -xyr_OpDep2(:,j,k) * xy_cosSZAInv(:,j) )
        end do
      end do

      !
      ! Eddington approximation
      !
      do k = 1, kmax
        do j = js, je
          xyz_Gam1(:,j,k) =  ( 7.0d0 - xyz_SSA2(:,j,k) * ( 4.0d0 + 3.0d0 * xyz_AF2(:,j,k) ) ) / 4.0d0
          xyz_Gam2(:,j,k) = -( 1.0d0 - xyz_SSA2(:,j,k) * ( 4.0d0 - 3.0d0 * xyz_AF2(:,j,k) ) ) / 4.0d0
          xyz_Gam3(:,j,k) =  ( 2.0d0 - 3.0d0 * xyz_AF2(:,j,k) * xy_cosSZA(:,j) )              / 4.0d0
          xyz_Gam4(:,j,k) = 1.0d0 - xyz_Gam3(:,j,k)
        end do
      end do


!!$    else if( sw .eq. 2 ) then
!!$
!!$         !
!!$         ! In infrared, delta-function adjustment is not done. 
!!$         !
!!$         af2  = af
!!$         ssa2 = ssa
!!$
!!$         do k = 1, km+1
!!$            do ij = ijs, ije
!!$               opdep2( ij, 1, k ) = opdep( ij, 1, k )
!!$            end do
!!$         end do
!!$
!!$         do k = 1, km+1
!!$            do ij = ijs, ije
!!$               exp_opdep2( ij, 1, k ) = 1.0d100
!!$            end do
!!$         end do
!!$
!!$         !
!!$         ! Hemispheric mean approximation
!!$         !
!!$         do k = 1, km
!!$            do ij = ijs, ije
!!$               gam1( ij, 1, k ) = 2.0d0 - ssa2 * ( 1.0d0 + af2 )
!!$               gam2( ij, 1, k ) = ssa2 * ( 1.0d0 - af2 )
!!$               gam3( ij, 1, k ) = 1.0d100
!!$               gam4( ij, 1, k ) = 1.0d100
!!$            end do
!!$         end do

      else
         stop 'Unexpected sw number in twostreamapp.f90'
      end if


      do k = 1, kmax
        do j = js, je
          xyz_DTau(:,j,k) = xyr_OpDep2(:,j,k-1) - xyr_OpDep2(:,j,k)
        end do
      end do
      !
      ! Avoiding singularity when dtau equal to zero 
      !
      do k = 1, kmax
        do j = js, je
          do i = 0, imax-1
            if( ( IDScheme .eq. IDSchemeLongWave ) .and. ( xyz_DTau(i,j,k) .lt. 1.0d-10 ) ) then
              xyz_DTau(i,j,k) = 0.0d0
            end if
          end do
        end do
      end do


      do k = 1, kmax
        do j = js, je
          xyz_Lambda(:,j,k) = sqrt( xyz_Gam1(:,j,k) * xyz_Gam1(:,j,k) - xyz_Gam2(:,j,k) * xyz_Gam2(:,j,k) )
          xyz_LSigma(:,j,k) = xyz_Gam2(:,j,k) / ( xyz_Gam1(:,j,k) + xyz_Lambda(:,j,k) )
        end do
      end do

      do k = 1, kmax
        do j = js, je
          xy_tmpval(:,j)       = exp( - xyz_Lambda(:,j,k) * xyz_DTau(:,j,k) )
          xyaz_smalle(:,j,1,k) = 1.0d0 + xyz_LSigma(:,j,k) * xy_tmpval(:,j)
          xyaz_smalle(:,j,2,k) = 1.0d0 - xyz_LSigma(:,j,k) * xy_tmpval(:,j)
          xyaz_smalle(:,j,3,k) = xyz_LSigma(:,j,k) + xy_tmpval(:,j)
          xyaz_smalle(:,j,4,k) = xyz_LSigma(:,j,k) - xy_tmpval(:,j)
        end do
      end do


      if( IDScheme .eq. IDSchemeShortWave ) then
        do k = 1, kmax
          do j = js, je
            xyz_cupb(:,j,k) = xyz_SSA2(:,j,k) * SolarFluxTOA * xyr_Trans2(:,j,k-1) * ( ( xyz_Gam1(:,j,k) - xy_cosSZAInv(:,j) ) * xyz_Gam3(:,j,k) + xyz_Gam2(:,j,k) * xyz_Gam4(:,j,k) ) / ( xyz_Lambda(:,j,k) * xyz_Lambda(:,j,k) - xy_cosSZAInvsq(:,j) )
            xyz_cupt(:,j,k) = xyz_SSA2(:,j,k) * SolarFluxTOA * xyr_Trans2(:,j,k  ) * ( ( xyz_Gam1(:,j,k) - xy_cosSZAInv(:,j) ) * xyz_Gam3(:,j,k) + xyz_Gam2(:,j,k) * xyz_Gam4(:,j,k) ) / ( xyz_Lambda(:,j,k) * xyz_Lambda(:,j,k) - xy_cosSZAInvsq(:,j) )
            xyz_cdob(:,j,k) = xyz_SSA2(:,j,k) * SolarFluxTOA * xyr_Trans2(:,j,k-1) * ( ( xyz_Gam1(:,j,k) + xy_cosSZAInv(:,j) ) * xyz_Gam4(:,j,k) + xyz_Gam2(:,j,k) * xyz_Gam3(:,j,k) ) / ( xyz_Lambda(:,j,k) * xyz_Lambda(:,j,k) - xy_cosSZAInvsq(:,j) )
            xyz_cdot(:,j,k) = xyz_SSA2(:,j,k) * SolarFluxTOA * xyr_Trans2(:,j,k  ) * ( ( xyz_Gam1(:,j,k) + xy_cosSZAInv(:,j) ) * xyz_Gam4(:,j,k) + xyz_Gam2(:,j,k) * xyz_Gam3(:,j,k) ) / ( xyz_Lambda(:,j,k) * xyz_Lambda(:,j,k) - xy_cosSZAInvsq(:,j) )
          end do
        end do

!!$      else if( sw .eq. 2 ) then
!!$
!!$         do k = 1, km
!!$            do ij = ijs, ije
!!$               !
!!$               ! Notice!
!!$               ! Avoiding singularity when dtau equal to zero. 
!!$               ! dtau occationally has much smaller value. 
!!$               ! When this occurs, b1 cannot be calculated correctly. 
!!$               !
!!$               if( dtau( ij, 1, k ) .ne. 0.0d0 ) then
!!$                  b0( ij, 1, k ) = pfinth( ij, 1, k )
!!$                  b1( ij, 1, k ) = ( pfinth( ij, 1, k+1 ) - pfinth( ij, 1, k ) ) / dtau( ij, 1, k )
!!$               else
!!$                  b0( ij, 1, k ) = 0.0d0
!!$                  b1( ij, 1, k ) = 0.0d0
!!$               end if
!!$
!!$            end do
!!$         end do
!!$
!!$         do k = 1, km
!!$            do ij = ijs, ije
!!$               mu1 = ( 1.0d0 - ssa2 ) / ( gam1( ij, 1, k ) - gam2( ij, 1, k ) )
!!$
!!$               cupb( ij, 1, k ) = pix2 * mu1 * ( b0( ij, 1, k ) &
!!$                    + b1( ij, 1, k ) &
!!$                    * ( dtau( ij, 1, k ) + 1.0d0 / ( gam1( ij, 1, k ) + gam2( ij, 1, k ) ) ) )
!!$               cupt( ij, 1, k ) = pix2 * mu1 * ( b0( ij, 1, k ) &
!!$                    + b1( ij, 1, k ) &
!!$                    * ( 0.0d0     + 1.0d0 / ( gam1( ij, 1, k ) + gam2( ij, 1, k ) ) ) )
!!$               cdob( ij, 1, k ) = pix2 * mu1 * ( b0( ij, 1, k ) &
!!$                    + b1( ij, 1, k ) &
!!$                    * ( dtau( ij, 1, k ) - 1.0d0 / ( gam1( ij, 1, k ) + gam2( ij, 1, k ) ) ) )
!!$               cdot( ij, 1, k ) = pix2 * mu1 * ( b0( ij, 1, k ) &
!!$                    + b1( ij, 1, k ) &
!!$                    * ( 0.0d0     - 1.0d0 / ( gam1( ij, 1, k ) + gam2( ij, 1, k ) ) ) )
!!$            end do
!!$         end do
!!$
      else
        stop 'Unexpected sw number in twostreamapp.f90'
      end if


      k = 1
      l = 1
      do j = js, je
        do i = 0, imax-1
          aa_TridiagMtx1( i+imax*(j-1)+1, l ) = 0.0d0
          aa_TridiagMtx2( i+imax*(j-1)+1, l ) = xyaz_smalle(i,j,2,k) - xy_SurfAlbedo(i,j) * xyaz_smalle(i,j,4,k)
          aa_TridiagMtx3( i+imax*(j-1)+1, l ) = xyaz_smalle(i,j,1,k) - xy_SurfAlbedo(i,j) * xyaz_smalle(i,j,3,k)
        end do
      end do
      if( IDScheme .eq. IDSchemeShortWave ) then
        do j = js, je
          do i = 0, imax-1
            aa_Vec( i+imax*(j-1)+1, l ) = - xyz_cupb(i,j,k) + xy_SurfAlbedo(i,j) * xyz_cdob(i,j,k) + xy_SurfAlbedo(i,j) * SolarFluxTOA * xyr_Trans2(i,j,0) * xy_cosSZA(i,j)
          end do
        end do
!!$      else if( sw .eq. 2 ) then
!!$         do ij = ijs, ije
!!$!            gemis = 1.0d0
!!$            gemis = emis( ij, 1 )
!!$            ea( (ij-ijs+1), l ) &
!!$                 = -cupb( ij, 1, km ) + ( 1.0d0 - gemis ) * cdob( ij, 1, km ) &
!!$                 + gemis * pi * pfinth( ij, 1, km+1 )
!!$         end do
      else
        stop 'Unexpected sw number in twostreamapp.f90'
      end if


!!$      do k = 1+1, kmax
!!$
!!$        do j = 1, jmax
!!$          do i = 0, imax-1
!!$
!!$            ! for even l
!!$            !
!!$            l = 2 * (k-1)
!!$            !
!!$            aa( i+imax*(j-1), l ) =                    &
!!$              &   smalle(i,j,3,k) * smalle(i,j,4,k-1)  &
!!$              & - smalle(i,j,1,k) * smalle(i,j,2,k-1)
!!$            ba( i+imax*(j-1), l ) =                    &
!!$              &   smalle(i,j,1,k) * smalle(i,j,1,k-1)  &
!!$              & - smalle(i,j,3,k) * smalle(i,j,3,k-1)
!!$            da( i+imax*(j-1), l ) =                    &
!!$              &   smalle(i,j,2,k) * smalle(i,j,3,k  )  &
!!$              & - smalle(i,j,1,k) * smalle(i,j,4,k  )
!!$            ea( i+imax*(j-1), l ) =                    &
!!$              &   smalle(i,j,3,k) * ( cupt(i,j,k-1) - cupb(i,j,k) ) &
!!$              & - smalle(i,j,1,k) * ( cdot(i,j,k-1) - cdob(i,j,k) )
!!$
!!$            ! for odd l
!!$            !
!!$            l = 2 * (k-1) + 1
!!$
!!$            aa( i+imax*(j-1), l ) = &
!!$              &   smalle(i,j,1,k-1) * smalle(i,j,4,k-1) &
!!$              & - smalle(i,j,2,k-1) * smalle(i,j,3,k-1)
!!$            ba( i+imax*(j-1), l ) = &
!!$              &   smalle(i,j,2,k  ) * smalle(i,j,2,k-1) &
!!$              & - smalle(i,j,4,k  ) * smalle(i,j,4,k-1)
!!$            da( i+imax*(j-1), l ) = &
!!$              &   smalle(i,j,1,k  ) * smalle(i,j,2,k-1) &
!!$              & - smalle(i,j,3,k  ) * smalle(i,j,4,k-1)
!!$            ea( i+imax*(j-1), l ) = &
!!$              &   smalle(i,j,2,k-1) * ( cupt(i,j,k-1) - cupb(i,j,k) ) &
!!$              & - smalle(i,j,4,k-1) * ( cdot(i,j,k-1) - cdob(i,j,k) )
!!$          end do
!!$        end do
!!$      end do


      do k = 1, kmax-1

        do j = js, je
          do i = 0, imax-1

            l = 2 * k     ! equation number
            !
            aa_TridiagMtx1( i+imax*(j-1)+1, l ) = xyaz_smalle(i,j,3,k+1) * xyaz_smalle(i,j,4,k  ) - xyaz_smalle(i,j,1,k+1) * xyaz_smalle(i,j,2,k  )
            aa_TridiagMtx2( i+imax*(j-1)+1, l ) = xyaz_smalle(i,j,1,k+1) * xyaz_smalle(i,j,1,k  ) - xyaz_smalle(i,j,3,k+1) * xyaz_smalle(i,j,3,k  )
            aa_TridiagMtx3( i+imax*(j-1)+1, l ) = xyaz_smalle(i,j,2,k+1) * xyaz_smalle(i,j,3,k+1) - xyaz_smalle(i,j,1,k+1) * xyaz_smalle(i,j,4,k+1)
            aa_Vec        ( i+imax*(j-1)+1, l ) = xyaz_smalle(i,j,3,k+1) * ( xyz_cupt(i,j,k  ) - xyz_cupb(i,j,k+1) ) - xyaz_smalle(i,j,1,k+1) * ( xyz_cdot(i,j,k  ) - xyz_cdob(i,j,k+1) )

            l = 2 * k + 1 ! equation number
            !
            aa_TridiagMtx1( i+imax*(j-1)+1, l ) = xyaz_smalle(i,j,1,k  ) * xyaz_smalle(i,j,4,k  ) - xyaz_smalle(i,j,2,k  ) * xyaz_smalle(i,j,3,k  )
            aa_TridiagMtx2( i+imax*(j-1)+1, l ) = xyaz_smalle(i,j,2,k+1) * xyaz_smalle(i,j,2,k  ) - xyaz_smalle(i,j,4,k+1) * xyaz_smalle(i,j,4,k  )
            aa_TridiagMtx3( i+imax*(j-1)+1, l ) = xyaz_smalle(i,j,1,k+1) * xyaz_smalle(i,j,2,k  ) - xyaz_smalle(i,j,3,k+1) * xyaz_smalle(i,j,4,k  )
            aa_Vec        ( i+imax*(j-1)+1, l ) = xyaz_smalle(i,j,2,k  ) * ( xyz_cupt(i,j,k  ) - xyz_cupb(i,j,k+1) ) - xyaz_smalle(i,j,4,k  ) * ( xyz_cdot(i,j,k  ) - xyz_cdob(i,j,k+1) )
          end do
        end do
      end do


      k = kmax
      l = 2 * kmax
      do j = js, je
        do i = 0, imax-1
          aa_TridiagMtx1( i+imax*(j-1)+1, l ) = -xyaz_smalle(i,j,2,k)
          aa_TridiagMtx2( i+imax*(j-1)+1, l ) =  xyaz_smalle(i,j,1,k)
          aa_TridiagMtx3( i+imax*(j-1)+1, l ) = 0.0d0
          aa_Vec        ( i+imax*(j-1)+1, l ) = -xyz_cdot(i,j,k) + 0.0d0
        end do
      end do

      ms = 0      + imax*(js-1)+1
      me = imax-1 + imax*(je-1)+1
      call tridiaginv( imax*jmax, 2*kmax, aa_TridiagMtx1, aa_TridiagMtx2, aa_TridiagMtx3, aa_Vec, ms, me )

      if( IDScheme .eq. IDSchemeShortWave ) then

        do k = 1, kmax
          do j = js, je
            do i = 0, imax-1
              xyr_UwFlux(i,j,k-1) = aa_Vec( i+imax*(j-1)+1, 2*k  ) * xyaz_smalle(i,j,1,k) + aa_Vec( i+imax*(j-1)+1, 2*k-1) * xyaz_smalle(i,j,2,k) + xyz_cupb(i,j,k)
              xyr_DwFlux(i,j,k-1) = aa_Vec( i+imax*(j-1)+1, 2*k  ) * xyaz_smalle(i,j,3,k) + aa_Vec( i+imax*(j-1)+1, 2*k-1) * xyaz_smalle(i,j,4,k) + xyz_cdob(i,j,k)
            end do
          end do
        end do

        k = kmax
        do j = js, je
          do i = 0, imax-1
            xyr_UwFlux(i,j,k) = aa_Vec( i+imax*(j-1)+1, 2*k  ) * xyaz_smalle(i,j,3,k) - aa_Vec( i+imax*(j-1)+1, 2*k-1) * xyaz_smalle(i,j,4,k) + xyz_cupt(i,j,k)
            xyr_DwFlux(i,j,k) = aa_Vec( i+imax*(j-1)+1, 2*k  ) * xyaz_smalle(i,j,1,k) - aa_Vec( i+imax*(j-1)+1, 2*k-1) * xyaz_smalle(i,j,2,k) + xyz_cdot(i,j,k)
          end do
        end do

        !
        ! Addition of Direct Solar Insident
        !
        do k = 0, kmax
          do j = js, je
            xyr_DwFlux(:,j,k) = xyr_DwFlux(:,j,k) + SolarFluxTOA * xyr_Trans2(:,j,k) * xy_cosSZA(:,j)
          end do
        end do

!!$      else if( sw .eq. 2 ) then
!!$
!!$         ! Source function technique described by Toon et al. [1989] 
!!$         ! is used to calculated infrared heating rate. 
!!$         !
!!$         do k = 1, km+1
!!$            do ij = ijs, ije
!!$               upfl( ij, 1, k ) = 0.0d0
!!$               dofl( ij, 1, k ) = 0.0d0
!!$            end do
!!$         end do
!!$
!!$         do l = 1, quadn
!!$            tmpcos = qucos( l )
!!$
!!$            k = 1
!!$            do ij = ijs, ije
!!$               intedo( ij, 1, k ) = 0.0d0
!!$            end do
!!$            do k = 1, km
!!$               do ij = ijs, ije
!!$                  tmpj = ( ea( (ij-ijs+1), k+k-1 ) + ea( (ij-ijs+1), k+k ) ) * lsigma( ij, 1, k ) &
!!$                       * ( lambda( ij, 1, k ) + 2.0d0 )
!!$                  tmpk = ( ea( (ij-ijs+1), k+k-1 ) - ea( (ij-ijs+1), k+k ) ) &
!!$                       * ( 2.0d0 - lambda( ij, 1, k ) )
!!$                  sig1 = pix2 * ( b0( ij, 1, k ) &
!!$                       - b1( ij, 1, k ) &
!!$                       * ( 1.0d0 / ( gam1( ij, 1, k ) + gam2( ij, 1, k ) ) - 0.5d0 ) )
!!$                  sig2 = pix2 * b1( ij, 1, k )
!!$                  intedo( ij, 1, k+1 ) = intedo( ij, 1, k ) &
!!$                       * exp( -dtau( ij, 1, k ) / tmpcos ) &
!!$                       + tmpj / ( lambda( ij, 1, k ) * tmpcos + 1.0d0 ) &
!!$                       * ( 1.0d0 &
!!$                       - exp( -dtau( ij, 1, k ) &
!!$                       * ( lambda( ij, 1, k ) + 1.0d0 / tmpcos ) ) ) &
!!$                       + tmpk / ( lambda( ij, 1, k ) * tmpcos - 1.0d0 ) &
!!$                       * ( exp( -dtau( ij, 1, k ) / tmpcos ) &
!!$                       -exp( -dtau( ij, 1, k ) * lambda( ij, 1, k ) ) ) &
!!$                       + sig1 * ( 1.0d0 - exp( -dtau( ij, 1, k ) / tmpcos ) ) &
!!$                       + sig2 * ( tmpcos * exp( -dtau( ij, 1, k ) / tmpcos ) &
!!$                       + dtau( ij, 1, k ) - tmpcos )
!!$               end do
!!$            end do
!!$
!!$            k = km+1
!!$            do ij = ijs, ije
!!$!               gemis = 1.0d0
!!$               gemis = emis( ij, 1 )
!!$               inteup( ij, 1, k ) = ( 1.0d0 - gemis ) * intedo( ij, 1, km+1 ) &
!!$                    + gemis * pix2 * pfinth( ij, 1, km+1 )
!!$            end do
!!$            do k = km, 1, -1
!!$               do ij = ijs, ije
!!$                  tmpg = ( ea( (ij-ijs+1), k+k-1 ) + ea( (ij-ijs+1), k+k ) ) &
!!$                       * ( 2.0d0 - lambda( ij, 1, k ) )
!!$                  tmph = ( ea( (ij-ijs+1), k+k-1 ) - ea( (ij-ijs+1), k+k ) ) * lsigma( ij, 1, k ) &
!!$                       * ( lambda( ij, 1, k ) + 2.0d0 )
!!$                  alp1 = pix2 * ( b0( ij, 1, k ) &
!!$                       + b1( ij, 1, k ) &
!!$                       * ( 1.0d0 / ( gam1( ij, 1, k ) + gam2( ij, 1, k ) ) - 0.5d0 ) )
!!$                  alp2 = pix2 * b1( ij, 1, k )
!!$                  inteup( ij, 1, k ) = inteup( ij, 1, k+1 ) &
!!$                       * exp( -dtau( ij, 1, k ) / tmpcos ) &
!!$                       + tmpg / ( lambda( ij, 1, k ) * tmpcos - 1.0d0 ) &
!!$                       * ( exp( -dtau( ij, 1, k ) / tmpcos ) &
!!$                       - exp( -dtau( ij, 1, k ) * lambda( ij, 1, k ) ) ) &
!!$                       + tmph / ( lambda( ij, 1, k ) * tmpcos + 1.0d0 ) &
!!$                       * ( 1.0d0 &
!!$                       -exp( -dtau( ij, 1, k ) &
!!$                       * ( lambda( ij, 1, k ) + 1.0d0 / tmpcos ) ) ) &
!!$                       + alp1 * ( 1.0d0 - exp( -dtau( ij, 1, k ) / tmpcos ) ) &
!!$                       + alp2 * ( tmpcos &
!!$                       - ( dtau( ij, 1, k ) + tmpcos ) &
!!$                       * exp( -dtau( ij, 1, k ) / tmpcos ) )
!!$               end do
!!$            end do
!!$
!!$            do k = 1, km+1
!!$               do ij = ijs, ije
!!$                  upfl( ij, 1, k ) = upfl( ij, 1, k ) &
!!$                       + inteup( ij, 1, k ) * qucos( l ) * qudcos( l )
!!$                  dofl( ij, 1, k ) = dofl( ij, 1, k ) &
!!$                       + intedo( ij, 1, k ) * qucos( l ) * qudcos( l )
!!$               end do
!!$            end do
!!$         end do
!!$
      else
        stop 'Unexpected sw number in twostreamapp.f90'
      end if


!!$      do ij = ijs, ije
!!$         gdf( ij, 1 ) = dofl( ij, 1, km+1 )
!!$      end do
!!$      if( sw == 1 ) then
!!$        ! visible
!!$        do ij = ijs, ije
!!$          gdf( ij, 1 ) = upfl( ij, 1, km+1 ) - dofl( ij, 1, km+1 )
!!$        end do
!!$      else
!!$        ! ir
!!$        do ij = ijs, ije
!!$          gdf( ij, 1 ) = dofl( ij, 1, km+1 )
!!$        end do
!!$      end if
!!$      do ij = ijs, ije
!!$        gdf( ij, 1 ) = upfl( ij, 1, km+1 ) - dofl( ij, 1, km+1 )
!!$      end do


!!$      do k = 1, km
!!$         do ij = ijs, ije
!!$            q( ij, 1, k ) &
!!$                 = ( ( upfl( ij, 1, k+1 ) - dofl( ij, 1, k+1 ) ) &
!!$                   - ( upfl( ij, 1, k   ) - dofl( ij, 1, k   ) ) ) &
!!$                   / ( gph ( ij, 1, k+1 ) - gph ( ij, 1, k   ) ) &
!!$                   * grav / cp
!!$         end do
!!$      end do


!!$      if( sw == 1 ) then
!!$        write( 6, * ) 'vis flux=', &
!!$!          -gdf((ijs+ije)/2,1), &
!!$!          -gdf((ijs+ije)/2,1) * ( 1.0d0 - albedo((ijs+ije)/2,1) ), &
!!$          -gdf((ijs+ije)/2,1) * albedo((ijs+ije)/2,1), &
!!$!          upfl((ijs+ije)/2,1,km+1)-dofl((ijs+ije)/2,1,km+1),
!!$          upfl((ijs+ije)/2,1,km+1), &
!!$          upfl((ijs+ije)/2,1,km+1) - gdf((ijs+ije)/2,1) * albedo((ijs+ije)/2,1)
!!$      else
!!$        write( 6, * ) 'ir  flux=', &
!!$          -gdf((ijs+ije)/2,1), &
!!$          upfl((ijs+ije)/2,1,km+1)-dofl((ijs+ije)/2,1,km+1), &
!!$          ' sig * Ts^4 = ', sboltz * gts((ijs+ije)/2,1)**4
!!$      end if

      do k = 0, kmax
        do j = js, je
          xyr_RadFlux(:,j,k) = xyr_UwFlux(:,j,k) - xyr_DwFlux(:,j,k)
        end do
      end do

      if( IDScheme .eq. IDSchemeShortWave ) then
        do j = js, je
          do i = 0, imax-1
            if( xy_cosSZA(i,j) <= 0.0d0 ) then
              do k = 0, kmax
                xyr_RadFlux(i,j,k) = 0.0d0
              end do
            end if
          end do
        end do
      end if


!!$      !
!!$      ! output variables
!!$      !
!!$      do ij = ijs, ije
!!$         goru( ij, 1 ) = upfl( ij, 1, 1 )
!!$         gord( ij, 1 ) = dofl( ij, 1, 1 )
!!$         gsru( ij, 1 ) = upfl( ij, 1, km+1 )
!!$         gsrd( ij, 1 ) = dofl( ij, 1, km+1 )
!!$      end do
!!$      if( sw .eq. 1 ) then
!!$         do ij = ijs, ije
!!$            if( cosz( ij, 1 ) .le. 0.0d0 ) then
!!$               goru( ij, 1 ) = 0.0d0
!!$               gord( ij, 1 ) = 0.0d0
!!$               gsru( ij, 1 ) = 0.0d0
!!$               gsrd( ij, 1 ) = 0.0d0
!!$            end if
!!$         end do
!!$      end if
!!$      do ij = ijs, ije
!!$         gor ( ij, 1 ) = goru( ij, 1 ) - gord( ij, 1 )
!!$         gsr ( ij, 1 ) = gsru( ij, 1 ) - gsrd( ij, 1 )
!!$      end do


    end subroutine RadiationTwoStreamAppCore
RadiationTwoStreamAppWrapper( QeRatio, xyz_SSA, xyz_AF, SolarFluxTOA, xy_SurfAlbedo, IDScheme, xy_InAngle, xyr_OptDepBase, xyr_RadFlux )
Subroutine :
QeRatio :real(DP), intent(in )
xyz_SSA( 0:imax-1, 1:jmax, 1:kmax ) :real(DP), intent(in )
xyz_AF( 0:imax-1, 1:jmax, 1:kmax ) :real(DP), intent(in )
SolarFluxTOA :real(DP), intent(in )
xy_SurfAlbedo( 0:imax-1, 1:jmax ) :real(DP), intent(in )
IDScheme :integer , intent(in )
xy_InAngle(0:imax-1, 1:jmax) :real(DP), intent(in )
: sec (入射角). sec (angle of incidence)
xyr_OptDepBase( 0:imax-1, 1:jmax, 0:kmax ) :real(DP), intent(in )
xyr_RadFlux( 0:imax-1, 1:jmax, 0:kmax ) :real(DP), intent(out)

[Source]

  subroutine RadiationTwoStreamAppWrapper( QeRatio, xyz_SSA, xyz_AF, SolarFluxTOA, xy_SurfAlbedo, IDScheme, xy_InAngle, xyr_OptDepBase, xyr_RadFlux )

    ! USE statements
    !

    ! 格子点設定
    ! Grid points settings
    !
    use gridset, only: imax, jmax, kmax    ! 鉛直層数.
                               ! Number of vertical level

    ! OpenMP
    !
    !$ use omp_lib


    real(DP), intent(in ) :: QeRatio
    real(DP), intent(in ) :: xyz_SSA       ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP), intent(in ) :: xyz_AF        ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP), intent(in ) :: SolarFluxTOA
    real(DP), intent(in ) :: xy_SurfAlbedo ( 0:imax-1, 1:jmax )
    integer , intent(in ) :: IDScheme
    real(DP), intent(in ) :: xy_InAngle    (0:imax-1, 1:jmax)
                              ! sec (入射角).
                              ! sec (angle of incidence)
    real(DP), intent(in ) :: xyr_OptDepBase( 0:imax-1, 1:jmax, 0:kmax )
    real(DP), intent(out) :: xyr_RadFlux   ( 0:imax-1, 1:jmax, 0:kmax )


    ! Local variables
    !
    integer :: js
    integer :: je

    integer :: nthreads
    integer, allocatable :: a_js(:)
    integer, allocatable :: a_je(:)

    integer :: n


    nthreads = 1
    !$ nthreads  = omp_get_max_threads()
!!$    !$ write( 6, * ) "Number of processors : ", omp_get_num_procs()
!!$    !$ write( 6, * ) "Number of threads    : ", nthreads

    allocate( a_js(0:nthreads-1) )
    allocate( a_je(0:nthreads-1) )

    do n = 0, nthreads-1

      if ( n == 0 ) then
        a_js(n) = 1
      else
        a_js(n) = a_je(n-1) + 1
      end if

      a_je(n) = a_js(n  ) + jmax / nthreads - 1
      if ( n + 1 <= mod( jmax, nthreads ) ) then
        a_je(n) = a_je(n) + 1
      end if

    end do


    !$OMP PARALLEL DEFAULT(PRIVATE) &
    !$OMP SHARED(nthreads,a_js,a_je, &
    !$OMP        QeRatio, xyz_SSA, xyz_AF, &
    !$OMP        SolarFluxTOA, &
    !$OMP        xy_SurfAlbedo, &
    !$OMP        IDScheme, &
    !$OMP        xy_InAngle, &
    !$OMP        xyr_OptDepBase, &
    !$OMP        xyr_RadFlux)

    !$OMP DO

    do n = 0, nthreads-1

      js = a_js(n)
      je = a_je(n)

      if ( js > je ) cycle

!!$      write( 6, * ) n, js, je

      call RadiationTwoStreamAppCore( QeRatio, xyz_SSA, xyz_AF, SolarFluxTOA, xy_SurfAlbedo, IDScheme, xy_InAngle, xyr_OptDepBase, xyr_RadFlux, js, je )

    end do


    !$OMP END DO
    !$OMP END PARALLEL


    deallocate( a_js )
    deallocate( a_je )



  end subroutine RadiationTwoStreamAppWrapper
module_name
Constant :
module_name = ‘radiation_two_stream_app‘ :character(*), parameter
: モジュールの名称. Module name
tridiag( mm, jmx, a, b, c, f, ms, me )
Subroutine :
mm :integer , intent(in )
jmx :integer , intent(in )
a( mm, jmx ) :real(DP), intent(in )
b( mm, jmx ) :real(DP), intent(in )
c( mm, jmx ) :real(DP), intent(in )
f( mm, jmx ) :real(DP), intent(inout)
ms :integer , intent(in )
me :integer , intent(in )

[Source]

    subroutine tridiag( mm, jmx, a, b, c, f, ms, me )

      integer , intent(in   ) :: mm, jmx
      real(DP), intent(in   ) :: a( mm, jmx ),b( mm, jmx )
      real(DP), intent(in   ) :: c( mm, jmx )
      real(DP), intent(inout) :: f( mm, jmx )
      integer , intent(in   ) :: ms, me


      ! Local variables
      !
      real(DP) :: q( mm, jmx ), p
      integer  :: j, m


      ! Forward elimination sweep
      !
      do m = ms, me
        q( m, 1 ) = - c( m, 1 ) / b( m, 1 )
        f( m, 1 ) =   f( m, 1 ) / b( m, 1 )
      end do

      do j = 2, jmx
        do m = ms, me
          p         = 1.0d0 / ( b( m, j ) + a( m, j ) * q( m, j-1 ) )
          q( m, j ) = - c( m, j ) * p
          f( m, j ) = ( f( m, j ) - a( m, j ) * f( m, j-1 ) ) * p
        end do
      end do

      ! Backward pass
      !
      do j = jmx - 1, 1, -1
        do m = ms, me
          f( m, j ) = f( m, j ) + q( m, j ) * f( m, j+1 )
        end do
      end do

    end subroutine tridiag
tridiag1( jmx, a, b, c, f )
Subroutine :
jmx :integer , intent(in )
a(jmx) :real(DP), intent(in )
b(jmx) :real(DP), intent(in )
c(jmx) :real(DP), intent(in )
f(jmx) :real(DP), intent(inout)

[Source]

    subroutine tridiag1( jmx, a, b, c, f )

      integer , intent(in   ) :: jmx
      real(DP), intent(in   ) :: a(jmx),b(jmx)
      real(DP), intent(in   ) :: c(jmx)
      real(DP), intent(inout) :: f(jmx)


      ! Local variables
      !
      real(DP) :: q(jmx), p
      integer  :: j


      ! Forward elimination sweep
      !
      q( 1 ) = - c( 1 ) / b( 1 )
      f( 1 ) =   f( 1 ) / b( 1 )

      do j = 2, jmx
        p      = 1.0d0 / ( b( j ) + a( j ) * q( j-1 ) )
        q( j ) = - c( j ) * p
        f( j ) = ( f( j ) - a( j ) * f( j-1 ) ) * p
      end do

      ! Backward pass
      !
      do j = jmx - 1, 1, -1
        f( j ) = f( j ) + q( j ) * f( j+1 )
      end do

    end subroutine tridiag1
tridiaginv( mm, jmx, a, b, c, f, ms, me )
Subroutine :
mm :integer , intent(in )
jmx :integer , intent(in )
a( mm, jmx ) :real(DP), intent(in )
b( mm, jmx ) :real(DP), intent(in )
c( mm, jmx ) :real(DP), intent(in )
f( mm, jmx ) :real(DP), intent(inout)
ms :integer , intent(in )
me :integer , intent(in )

[Source]

    subroutine tridiaginv( mm, jmx, a, b, c, f, ms, me )

      integer , intent(in   ) :: mm, jmx
      real(DP), intent(in   ) :: a( mm, jmx ),b( mm, jmx )
      real(DP), intent(in   ) :: c( mm, jmx )
      real(DP), intent(inout) :: f( mm, jmx )
      integer , intent(in   ) :: ms, me


      ! Local variables
      !
      real(DP) :: a2( mm, jmx )
      real(DP) :: b2( mm, jmx )
      real(DP) :: c2( mm, jmx )
      real(DP) :: f2( mm, jmx )
      integer  :: j
      integer  :: m

      do j = 1, jmx
        do m = ms, me
          a2(m,j) = c(m,jmx-j+1)
          b2(m,j) = b(m,jmx-j+1)
          c2(m,j) = a(m,jmx-j+1)
          f2(m,j) = f(m,jmx-j+1)
        end do
      end do

      call tridiag( mm, jmx, a2, b2, c2, f2, ms, me )

      do j = 1, jmx
        do m = ms, me
          f(m,j) = f2(m,jmx-j+1)
        end do
      end do

    end subroutine tridiaginv
version
Constant :
version = ’$Name: dcpam5-20100424 $’ // ’$Id: radiation_two_stream_app.f90,v 1.2 2010-03-24 12:10:47 yot Exp $’ :character(*), parameter
: モジュールのバージョン Module version