기상 기후 강의 노트

 ! File: solve_heat_diffusion_v1.f90
! Purpose: OOP implementation in Fortran of the ADE method of Barakat & Clark
!          (1966), for solving the time-dependent heat-diffusion equation in 2D,
 !          with a linearly-varying temperature profile on the domain-boundaries.
 !NOTES: - (for gfortran users): gfortran-4.8 still had no support for
 !         'final'-methods, and an error will be raised by that version (or
 !         earlier) of that compiler; HOWEVER, the final-procedure can be safely
 !         removed in our case, since it is included for demonstration-purposes
  !         only (the arrays which it de-allocates are freed anyway by the runtime
  !         system, when the 'Solver'-instance goes out of scope). Comment
  !         corresponding line below, to disable use of this feature.
 !       - The file "config_file_formatted.in", containing the parameters for the
 !         solver, needs to be located in working directory.
 !       - Lines containing only comments with dots can be ignored (they just
 !         help with including code in TeX-files).
 
 17 module NumericKinds
 18   implicit none
 
    ! KINDs for different types of REALs
 21   integer, parameter :: &
        R_SP = selected_real_kind(  6,   37 ), &
        R_DP = selected_real_kind( 15,  307 ), &
       R_QP = selected_real_kind( 33, 4931 )
   ! Alias for precision that we use in the program (change this to any of the
   ! values 'R_SP', 'R_DP', or 'R_QP', to switch to another precision globally).
 27   integer, parameter :: RK = R_DP ! if changing this, also change RK_FMT
 
    ! KINDs for different types of INTEGERs
 30   integer, parameter :: &
         I1B = selected_int_kind(2), & ! max = 127
         I2B = selected_int_kind(4), & ! max ~ 3.28x10^4
         I3B = selected_int_kind(9), & ! max ~ 2.15x10^9
         I4B = selected_int_kind(18)   ! max ~ 9.22x10^18
    ! Alias for integer-precision (analogue role to RK above).
 36   integer, parameter :: IK = I3B

    ! Edit-descriptors for real-values
 39   character(len=*), parameter :: R_SP_FMT = "f0.6", &
         R_DP_FMT = "f0.15", R_QP_FMT = "f0.33"
    ! Alias for output-precision to use in the program (keep this in sync with RK)
 42   character(len=*), parameter :: RK_FMT = R_DP_FMT
 43 end module NumericKinds

 45 module Config_class
 46   use NumericKinds
 47   implicit none
 48   private

 50   type, public :: Config
 51      real(RK) :: mDiffusivity = 1.15E-6_RK, & ! sandstone
            ! NOTE: "physical" units here (Celsius)
	      mTempA = 100._RK, &
            mTempB =  75._RK, &
            mTempC =  50._RK, &
            mTempD =  25._RK, &
			          mSideLength = 30._RK
 58      integer(IK) :: mNx = 200 ! # of points for square side-length
 59   end type Config

    ! Generic IFACE for user-defined CTOR
 62   interface Config
 63      module procedure createConfig
 64   end interface Config

 66 contains
 67   type(Config) function createConfig( cfgFilePath )
 68     character(len=*), intent(in) :: cfgFilePath
 69     integer :: cfgFileID
 70     open( newunit=cfgFileID, file=trim(cfgFilePath), status='old', action='read' )
      ! .....................................................
 72     read(cfgFileID,*) createConfig%mTempA
 73     read(cfgFileID,*) createConfig%mTempB
 74     read(cfgFileID,*) createConfig%mTempC
 75     read(cfgFileID,*) createConfig%mTempD
 76     read(cfgFileID,*) createConfig%mNx
 77     read(cfgFileID,*) createConfig%mDiffusivity
 78     read(cfgFileID,*) createConfig%mSideLength
 79     close(cfgFileID)
 80   end function createConfig
 81 end module Config_class
 
 83 module Solver_class
 84   use NumericKinds
 85   use Config_class
 86   implicit none
 87   private

 89   type, public :: Solver
 90      private ! Hide internal-data from users.
 91      type(Config) :: mConfig
 92      real(RK) :: mNt, & ! # of iterations to simulate a characteristic time
            mDx, mDt, mA, mB ! Configuration-dependent factors.
 94      real(RK), allocatable, dimension(:,:) :: mU, mV ! main work-arrays
 95      integer(IK) :: mNumItersMax, mCurrIter = 0
 96    contains
 97      private ! By default, hide methods (and expose as needed).
 98      procedure, public :: init
 99      procedure, public :: run
100      procedure, public :: writeAscii
101      procedure, public :: getTemp
 ! Internal methods (users don't need to know about these).
103      procedure :: advanceU
104      procedure :: advanceV
    !final :: cleanup ! NOTE: may need to comment-out for gfortran!
106   end type Solver

108	contains
 subroutine init( this, cfgFilePath, simTime ) ! initialization subroutine
110     class(Solver), intent(inout) :: this
111     character(len=*), intent(in) :: cfgFilePath
112     real(RK), intent(in) :: simTime
  ! .....................................................
114     integer(IK) :: nX, i, j
115     real(RK) :: lambda

117     this%mConfig = Config( cfgFilePath ) ! call component CTOR
118     nX = this%mConfig%mNx ! for making code below more compact
     ! conservative choice for N_t, to resolve transients -- see Barakat & Clark (1966)
120     this%mNt = nX**2
     ! evaluate derived parameter 'mLambda' in Solver, based on configuration
122     this%mNumItersMax = int( simTime*this%mNt )
123     this%mDx = 1. / nX
124     this%mDt = 1. / this%mNt
125     lambda = (2.*this%mDt) / (this%mDx**2)
126     this%mA = (1.-lambda)/(1.+lambda)
127     this%mB = lambda/(2.*(1.+lambda))
     ! allocate memory for internal arrays
129     allocate( this%mU(0:nX, 0:nX), this%mV(0:nX, 0:nX) )
     ! initialize mU-field:
     ! - set initial temperature everywhere...
132     this%mU = 1.
     ! - BUT re-write @ boundaries, for correct BCs
     ! -- North
135     this%mU(:, nX) = [ (1./3.*(i/real(nX, RK))+2./3., i=0,nX) ]
     ! -- West
137     this%mU(0, :) = [ (1./3.*(j/real(nX, RK))+1./3., j=0,nX) ]
     ! -- South
139     this%mU(:, 0) = [ (-1./3.*(i/real(nX, RK))+1./3., i=0,nX) ]
     ! -- East
141     this%mU(nX, :) = [ (j/real(nX, RK), j=0,nX) ]
     ! initialize mV-field (from mU-field)
143     this%mV = this%mU
144   end subroutine init

146   real(RK) function getTemp( this, i, j ) ! GETter for temperature
147     class(Solver), intent(in) :: this
148     integer(IK), intent(in) :: i, j
149     getTemp = 0.5*( this%mU(i,j) + this%mV(i,j) )
150	end function getTemp

	
	 subroutine run( this ) ! method for time-marching
153     class(Solver), intent(inout) :: this
154     integer(IK) :: k ! dummy index (time-marching)

156     do k=1, this%mNumItersMax ! MAIN loop
       ! simple progress-monitor
158        if( mod(k-1, (this%mNumItersMax-1)/10) == 0 ) then
159           write(*, '(i5,a)') nint((k*100.0)/this%mNumItersMax), "%"
160        end if
        ! defer work to private methods
162        call this%advanceU()
163        call this%advanceV()
164        this%mCurrIter = this%mCurrIter + 1 ! tracking time step
165     end do
166   end subroutine run

168   subroutine advanceU( this )
169     class(Solver), intent(inout) :: this
170     integer(IK) :: i, j ! local variables
 ! actual update for 'mU'-field (NE-ward)
172     do j=1, this%mConfig%mNx-1   ! do NOT update
173        do i=1, this%mConfig%mNx-1 ! boundaries
174           this%mU(i,j) = this%mA*this%mU(i,j) + this%mB*( &
                this%mU(i-1,j) + this%mU(i+1,j) + this%mU(i,j-1) + this%mU(i,j+1) )
176        end do
177     end do
178   end subroutine advanceU

180   subroutine advanceV( this ) ! similar to 'advanceU'
181     class(Solver), intent(inout) :: this
    ! .....................................................
183     integer(IK) :: i, j ! local variables
    ! actual update for 'mV'-field (SW-ward)
185     do j=this%mConfig%mNx-1, 1, -1   ! do NOT update
186        do i=this%mConfig%mNx-1, 1, -1 ! boundaries
187           this%mV(i,j) = this%mA*this%mV(i,j) + this%mB*( &
               this%mV(i-1,j) + this%mV(i+1,j) + this%mV(i,j-1) + this%mV(i,j+1) )
189        end do
190     end do
191	end subroutine advanceV

	 ! method for producing a ASCII output file
194   subroutine writeAscii( this, outFilePath )
195     class(Solver), intent(in) :: this
196     character(len=*), intent(in) :: outFilePath
    ! .....................................................
198     integer(IK) :: x, y, outFileID ! local variables

200     open( newunit=outFileID, file=trim(outFilePath), status='replace', action='write' )
201     write(outFileID, '(a)') &
          "# output file for program solve_heat_diffusion_v1.f90"
203     write(outFileID, '(a,2x,a)') '"s"', "# time unit"
204     write(outFileID, '(f0.8,2x,a)') &
         (this%mCurrIter*this%mConfig%mSideLength**2)/ &
          (this%mConfig%mDiffusivity*this%mNt), &
          "# current time"
208     write(outFileID, '(a,2x,a)') '"m"', "# X unit"
209     write(outFileID, '(i0,2x,a)') this%mConfig%mNx, "# Nx"
210     write(outFileID, '(a,2x,a)') '"m"', "# Y unit"
211     write(outFileID, '(i0,2x,a)') this%mConfig%mNx, "# Ny"
212     write(outFileID, '(a,2x,a)') '"degree~C"', "# temperature unit"
     ! X-axis
214     do x=0, this%mConfig%mNx
215        write(outFileID, '(f0.8,2x)', advance='no') this%mDx*this%mConfig%mSideLength*x
216     end do
217     write(outFileID, '(a)') "# XVals"
    ! Y-axis
219     do y=0, this%mConfig%mNx
220        write(outFileID, '(f0.8,2x)', advance='no') this%mDx*this%mConfig%mSideLength*y
221     end do
222     write(outFileID, '(a)') "# YVals"
    ! simulation results
224     write(outFileID, '(a)') "# from next line to end: simulated temperature"
225     do y=0, this%mConfig%mNx
226        do x=0, this%mConfig%mNx
227           write(outFileID, '(f0.8,2x)', advance='no') &
               this%mConfig%mTempD+this%getTemp(x,y)*(this%mConfig%mTempA-this%mConfig%mTempD)
229        end do
230        write(outFileID, *) ! newline to separate rows for R visualization script
231     end do
232     close(outFileID)
233   end subroutine writeAscii

  ! destructor method
236   subroutine cleanup( this )
     ! 'class' -> 'type' (dummy-arg cannot be polymorphic for final procedures)
238     type(Solver), intent(inout) :: this
    ! in this version, we only deallocate memory
240     deallocate( this%mU, this%mV )
241   end subroutine cleanup
242	end module Solver_class

	
	 program solve_heat_diffusion_v1
   use NumericKinds
   use Solver_class
247   implicit none

 type(Solver) :: square
   real(RK) :: simTime = 0.1 ! no. of characteristic time-intervals to simulate

252   character(len=200) :: configFile = "config_file_formatted.in", &
        outputFile = "simulation_final_temp_field.dat"

255   call square%init( configFile, simTime ) ! call Initializer
256   call square%run()

258	  call square%writeAscii( outputFile )
	  
	  pause
259 end program solve_heat_diffusion_v1

! 온도변환 인터페이스 사용 버전 p222
PROGRAM Temperature_Conversion_7
    implicit none
    
    INTERFACE
        Function Celsius_to_Fahr(Temp)        
        real :: Celsius_to_Fahr
        real, intent(in) :: Temp                
        END Function Celsius_to_Fahr
    END INTERFACE
    
    real :: fahrenheit, celsius
    character(1) :: response
        
    DO
        ! Get a Celsius temperature
        write (*, '(1x, A)', ADVANCE = "NO") "Enter a Celsius temperature:"
        read *, Celsius
    
        ! Use the module function Fahr_to_Celsius to convert it to Celsius       
        Fahrenheit = Celsius_to_Fahr(celsius)
        
        ! Output the result
        print '(1x, 2(F6.2, A))', celsius, & 
            " in Celsius is equivalent to ", fahrenheit, " in Fahrenheit"
        
        ! Check if more temperautre ar to ber converted
        write (*, '(/ 1x, A)', ADVANCE = "NO") &
               "More temperatures to convert (Y or N)?"
        read *, response
                
        IF (response /= "Y") EXIT
    END DO
pause
    END PROGRAM Temperature_Conversion_7

    
Function Celsius_to_Fahr(Temp)
    implicit none
    real :: Celsius_to_Fahr
    real, intent(in) :: Temp        
    Celsius_to_Fahr = (Temp - 32.0) /1.8
End Function Celsius_to_Fahr

https://aeir.tistory.com/entry/%ED%8F%AC%ED%8A%B8%EB%9E%80-%EA%B0%95%EC%A2%8C-%EC%84%9C%EB%B8%8C%EB%A3%A8%ED%8B%B4-%EB%B6%80%ED%94%84%EB%A1%9C%EA%B7%B8%EB%9E%A8?category=940076 

1. 오염지수 구하기 - IF 문

PROGRAM pollution
    IMPLICIT NONE
    
   INTEGER :: level_1, level_2, level_3, index
   INTEGER, PARAMETER :: cutoff = 50
   
   ! get the 3 pollution readings
   PRINT *, "Enter 3 pollution readings (unit: ppm) :"
   READ *, level_1, level_2, level_3
   
   index = (level_1 + level_2 + level_3)/3
   
   IF (index < cutoff) THEN 
       PRINT *, "SAFE"
       ELSE
       PRINT *, "HAZARDOUS"
   END IF
   
pause
    end program pollution

2. 오염지수 구하기 - IF-ELSE 문

program pollution
    implicit none
    
   integer :: level_1, level_2, level_3, index
   integer, parameter :: lowcutoff = 25, highcutoff = 50
   
   ! get the 3 pollution readings
   print *, "Enter 3 pollution readings (unit: ppm) :"
   read *, level_1, level_2, level_3
   
   ! Calculate the pollution index 
   index = (level_1 + level_2 + level_3)/3
   
   ! classify the pollution index and display air-quality conditions     
   if (index < lowcutoff) then 
       print *, "GOOD"       
   else if (index < highcutoff) then 
           print *, "FAIR"
   else 
           print *, "POOR"
   end if
       
   
   
pause
    end program pollution

https://aeir.tistory.com/entry/%ED%8F%AC%ED%8A%B8%EB%9E%80-%EB%85%BC%EB%A6%AC%EC%8B%9D?category=940076 

https://aeir.tistory.com/entry/%ED%8F%AC%ED%8A%B8%EB%9E%80-IF%EC%99%80-IF-ELSE-%EA%B5%AC%EB%AC%B8?category=940076 

과제

화씨와 섭씨를 구분해 변환하는 하나의 코드 작성

1. 125^1/3 승수 구하기

program ex
    implicit none
    
    real :: x, y1, y2, y3, y4,  y5
    x = 125    
    
    y1 = x ** (1.0/3.0)
    y2 = x ** (1./3.)
    y3 = x ** (1/3.0)
    y4 = x ** (1.0/3)
    y5 = x ** (1/3)
    print *, '12345678901234567890123456789012'
    print *, 'Remind 1.0 x 3.0 = ', 1.0 * 3.0
    print *, 'Remind   1 x   3 = ', 1 * 3
    print *
    
    print *, y1
    print *, y2
    print *, y3
    print *, y4
    print *, y5
    
    print *
pause
end program ex

결과해설

사칙 연산에서 아래와 같은 사실을 알 수 있다. 

실수 ÷ 실수 = 실수

실수 ÷ 정수 = 실수

정수 ÷ 실수 = 실수

정수 ÷ 정수 = 정수

y5의 경우 1/3 = 0.333 인데, 정수/정수=정수형이 되므로, 소수점 이하자리는 버림되고 1/3은 0이 된다. 따라서, 결과값이 1이다. 1은 정수가 아니라 1.00000으로 실수형으로 표현되는 이유는 y5가 실수형으로 선언되었기 때문이다. 

3. e³ 구하기

입력설계

키보드로 입력값을 받아들인다. 

내장함수 EXP( )를 사용하여 지수값을 구한다. 

의사코드

   read x

   p = EXP(x)

   print p

프로그램

program ex
    implicit none
    
    real :: x, p
    
    print *, 'Enter value :'
    read *, x
    
    p = exp(x)
    print *, 'Exponential value is', p
    
    print *
pause
end program ex

결과 해석

원하는 값을 입력받아 프로그램을 실행하고자 할 때, read 문을 사용한다. 

4. log_e2.7 와log₁₀2.7구하기

처리조건

자연로그 값을 구하기 위해 log( )와 alog( )함수를 사용해서 비교한다. 

상용로그 값을 구하기 위해 log10( )함수를 사용한다. 

프로그램

program ex
    implicit none
    
    real :: x, a, b, c
    
    print *, 'Enter value :'
    read *, x
    
    a = log(x)
    b = alog(x)
    c = log10(x)
    print *, ' LOG(x) is ', a
    print *, 'ALOG(x) is ', b
    print *, 'LOG10(x) is ', c
    print *
pause
    end program ex

5. 삼각함수 값 구하기

처리조건

  삼각함수 입력값은 radian 이다. degree를 radian 으로 변환해야 한다. 

program ex
    implicit none
    
    real :: degree, radian
    real :: x, L, M, N, P
    
    print *, 'Enter the angle in degree :'
    read *, degree
    
    radian = 3.141592/180.0
    x = radian * degree
    L = sin(x)
    M = cos(x)
    N = tan(x)
    P = sin(x)/cos(x)
        
    print *, "SIN(X) is ", L
    print *, "COS(X) is ", M
    print *, "TAN(X) is ", N
    print *, "SIN(X)/COS(X) is ", P
        
    
    print *
pause
    end program ex

결과 해설

마지막에 TAN(x)와 SIN(x)/COS(x) 값을 비교해 보면, 맨 마지막 소수자리 숫자가 다르다. 이는 truncation 문제로 인해 발생하는 것이다. 

과제

아래 열거된 내장함수를 모두 사용하여 코딩하여 결과를 발표한다.

처리조건:

입력값은 화면에서 임의로 입력받도록 한다. 

출력은 소수점 3째자리까지 출력되도록 한다. 

발표내용:

코딩 설계 및 구성을 설명한다. 

코딩을 line-by-line 으로 설명한다. 

여러가지 내장함수