Subversion Repositories lagranto.ecmwf

Rev

Go to most recent revision | Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
3 michaesp 1
module evaluate
2
 
3
! The user can assign values to parameters that can be used in expressions with 
4
! the subroutine defparam. The calling syntax is
5
!
6
!     call defparam(symbol,value) or call defparam(symbol,expr)
7
!
8
! where symbol is the desired parameter name; value is a real, integer, or 
9
! complex variable (single or double precision); and expr is a string
10
! containing an expression to be evaluated. The value obtained by evaluating the
11
! expression expr is associated with the parameter symbol. Parameter names must
12
! begin with a letter (a-z, A-Z) and must not be longer than 24 characters.
13
! Parameter names are not case dependent.
14
!
15
! An expression can be evaluated with the subroutine evalexpr. The calling
16
! syntax is 
17
!
18
!          call evalexpr(expr,value)
19
!
20
! where expr is a string containing the expression to be evaluated; value is the
21
! result (single or double precision real, complex or integer). The 
22
! expression can contain the arithmetic operations +, -, *, /, or ^ as well as
23
! the functions sin, cos, tan, log, ln, abs, exp, sqrt, real, imag, conjg, and 
24
! ang (the function ang calculates the phase angle of its complex argument). The
25
! expression can also contain numerical values and previously defined parameters
26
! Grouping by nested levels of parentheses is also allowed. The parameters pi
27
! and i (imaginary unit) are predefined. Complex numbers can be entered as a+i*b 
28
! if the parameter i has not been redefined by the user. Complex numbers can also
29
! be entered using complex(a,b).
30
! Example expression: 
31
!
32
!          conjg(((cos(x) + sqrt(a+i*b))^2+complex(ln(1.6e-4),20))/2)
33
!
34
! An equation of the form <symbol> = <expression> can be evaluated using the 
35
! subroutine evaleqn. The calling syntax is
36
!
37
!          call evaleqn(eqn)       
38
!
39
! where eqn is a string containing the equation. The right-hand-side of the
40
! equation is evaluated and assigned to the symbol given by the left-hand-side.
41
!
42
! The value assigned to a symbol can be retrieved using the subroutine getparam.
43
! The calling syntax is
44
!
45
!          call getparam(sym,value)
46
!
47
! where sym is a symbol string; value is a numeric variable (any of the six 
48
! standard types).
49
!
50
! The symbols and their values in the symbol table can be listed using the
51
! subroutine listvar. The variable ierr is always available following a call
52
! to any of the above subroutines and is zero if there were no errors. The
53
! possible nonzero values for ierr are
54
!
55
! 1       Expression empty
56
! 2       Parentheses don't match
57
! 3       Number string does not correspond to a valid number
58
! 4       Undefined symbol
59
! 5       Less than two operands for binary operation
60
! 6       No operand for unary plus or minus operators
61
! 7       No argument(s) for function
62
! 8       Zero or negative real argument for logarithm
63
! 9       Negative real argument for square root
64
! 10      Division by zero
65
! 11      Improper symbol format
66
! 12      Missing operator
67
! 13      Undefined function
68
! 14      Argument of tangent function a multiple of pi/2
69
!
70
use precision
71
use strings
72
 
73
save
74
private
75
public :: valuep,evalexpr,defparam,evaleqn,getparam,listvar,ierr
76
 
77
type item                         
78
  character(len=24):: char
79
  character :: type
80
end type item
81
 
82
type param                        
83
  character (len=24):: symbol
84
  complex(kc8):: value
85
end type param
86
 
87
interface defparam                
88
  module procedure strdef       ! value given by expression      
89
  module procedure valdef_dc    ! Double precision complex value
90
  module procedure valdef_sc    ! Single precision complex value
91
  module procedure valdef_dr    ! Double precision real value
92
  module procedure valdef_sr    ! Single precision real value
93
  module procedure valdef_di    ! Double precision integer value
94
  module procedure valdef_si    ! Single precision integer value
95
end interface
96
 
97
interface evalexpr
98
  module procedure evalexpr_dc  ! Double precision complex result
99
  module procedure evalexpr_sc  ! Single precision complex result
100
  module procedure evalexpr_dr  ! Double precision real result
101
  module procedure evalexpr_sr  ! Single precision real result
102
  module procedure evalexpr_di  ! Double precision integer result
103
  module procedure evalexpr_si  ! Single precision integer result
104
end interface
105
 
106
interface getparam
107
  module procedure getparam_dc  ! Double precision complex result
108
  module procedure getparam_sc  ! Single precision complex result
109
  module procedure getparam_dr  ! Double precision real result
110
  module procedure getparam_sr  ! Single precision real result
111
  module procedure getparam_di  ! Double precision integer result
112
  module procedure getparam_si  ! Single precision integer result
113
end interface
114
 
115
integer,parameter :: numtok=100  ! Maximum number of tokens 
116
type(param) :: params(100)       ! Symbol table
117
integer :: nparams=0,itop,ibin
118
complex(kc8) :: valstack(numtok) ! Stack used in evaluation of expression
119
type(item):: opstack(numtok)     ! Operator stack used in conversion to postfix
120
integer :: ierr                  ! Error flag
121
 
122
!**********************************************************************
123
 
124
contains
125
 
126
!**********************************************************************
127
 
128
 
129
SUBROUTINE EVALEXPR_DC(expr,val)    ! Evaluate expression expr for
130
                                    ! val double precision complex
131
 
132
character (len=*),intent(in) :: expr
133
complex(kc8) :: val
134
character (len=len(expr)+1) :: tempstr
135
character :: cop
136
integer :: isp(numtok)          ! On stack priority of operators in opstack
137
integer :: lstr
138
complex(kc8) :: cval,oper1,oper2
139
real(kr8) :: valr,vali
140
type(item):: token(numtok)      ! List of tokens ( a token is an operator or 
141
                                ! operand) in postfix order
142
type(item) :: x,junk,tok
143
 
144
ierr=0
145
token(1:)%char=' '
146
 
147
if(nparams == 0) then                  ! Initialize symbol table
148
  params(1)%symbol='PI'
149
  params(1)%value=(3.14159265358979_kr8,0.0_kr8)
150
  params(2)%symbol='I'
151
  params(2)%value=(0.0_kr8,1.0_kr8)
152
  nparams=2
153
end if
154
 
155
if(len_trim(expr) == 0) then           ! Expression empty
156
  ierr=1
157
  write(*,*) 'Error: expression being evaluated is empty'
158
  return
159
end if
160
 
161
tempstr=adjustl(expr)
162
call removesp(tempstr)   ! Removes spaces, tabs, and control characters
163
 
164
! ****************************************************************************
165
! STEP 1:  Convert string to token array. Each token is either an operator or
166
!          an operand. Token array will be in postfix (reverse Polish) order.
167
!*****************************************************************************
168
 
169
ntok=0
170
ibin=0
171
itop=0
172
do
173
  lstr=len_trim(tempstr)
174
  call get_next_token(tempstr(1:lstr),tok,icp,insp)
175
  select case(tok%type)
176
  case('S')
177
    ntok=ntok+1
178
    token(ntok)=tok
179
  case('E')
180
    do 
181
      if(itop < 1)exit
182
      call popop(x)        ! Output remaining operators on stack
183
      ntok=ntok+1
184
      token(ntok)=x
185
    end do
186
    ntok=ntok+1
187
    token(ntok)=tok
188
    exit
189
  case('R')  ! Token is right parenenthesis
190
    do 
191
      if(opstack(itop)%type == 'L') exit  ! Output operators on stack down
192
      call popop(x)                       ! to left parenthesis
193
      ntok=ntok+1
194
      token(ntok)=x
195
    end do                             
196
    call popop(junk)                      ! Remove left parenthesis from stack
197
    if(opstack(itop)%type == 'F') then    ! Output function name if present
198
      call popop(x)
199
      ntok=ntok+1
200
      token(ntok)=x
201
    end if                               
202
  case('D')  ! Token is comma
203
    do 
204
      if(opstack(itop)%type == 'L') exit  ! Output operators on stack down
205
      call popop(x)                       ! to left parenthesis
206
      ntok=ntok+1
207
      token(ntok)=x
208
    end do                              
209
  case('U','B','L','F') ! Token is operator, left parenthesis or function name
210
    do 
211
      if(isp(itop) < icp) exit            ! Output operators on stack having
212
      call popop(x)                       ! an instack priority that is
213
      ntok=ntok+1                         ! greater than or equal to the
214
      token(ntok)=x                       ! priority of the incoming operator  
215
    end do                                
216
    call pushop(tok)     ! Put incoming operator on stack                       
217
    isp(itop)=insp
218
  end select
219
end do
220
 
221
isum=0                                 ! Error check for matching parentheses
222
do i=1,ntok
223
  if(token(i)%type == 'L' ) isum=isum+1
224
  if(token(i)%type == 'R' ) isum=isum-1
225
end do
226
if(isum /= 0) then
227
  ierr=2
228
  write(*,*) 'Error in the evaluation of the expression ',trim(expr)
229
  write(*,*) "Parentheses don't match"
230
  write(*,*)
231
  return
232
end if
233
 
234
 
235
!*****************************************************************************
236
! STEP 2: Evaluate token string in postfix order
237
!*****************************************************************************
238
 
239
itop=0
240
do i=1,ntok
241
  x=token(i)
242
  select case(x%type)
243
  case('E')  ! Token is end token
244
    if(itop>1) then                
245
      ierr=12
246
      write(*,*) 'Error: missing operator in expression ',trim(expr)
247
      write(*,*)
248
      return
249
    end if
250
    call popval(val)               ! Final result left on stack of values
251
    exit
252
  case('S')  ! Token is operand
253
    call valuep(x%char,cval)       ! Evaluate operand
254
    if(ierr/=0) return
255
    call pushval(cval)             ! Put value of operand on stack
256
  case('B')  ! Token is a binary operator
257
    if(itop < 2) then
258
      ierr=5
259
      write(*,*) 'Error in evaluation of expression ',trim(expr)
260
      write(*,*) 'Less than two operands for binary operator  '&
261
                 ,trim(x%char)
262
      write(*,*)
263
      return
264
    end if                         
265
    call popval(oper1)             ! Pull off top two values from stack
266
    call popval(oper2)
267
    select case(trim(x%char))      ! Perform operation on values
268
    case('^')
269
      cval=oper2**oper1
270
    case('*')
271
      cval=oper2*oper1
272
    case('/')
273
      if(oper1 == (0._kr8,0._kr8)) then
274
        ierr=10
275
        write(*,*) 'Error in expression ',trim(expr)
276
        write(*,*) 'Division by zero'
277
        write(*,*)
278
        return
279
      end if
280
      cval=oper2/oper1
281
    case('+')
282
      cval=oper2+oper1
283
    case('-')
284
      cval=oper2-oper1
285
    end select
286
    call pushval(cval)             ! Put result back on stack
287
  case('U')  ! Token is unary operator
288
    if(itop == 0) then
289
      ierr=6
290
      write(*,*) 'Error in expression ',trim(expr)
291
      write(*,*) 'No operand for unary operator ',trim(x%char)
292
      write(*,*)
293
      return
294
    else
295
      call popval(oper1)           ! Pull top value off stack
296
    end if
297
    select case(trim(x%char))      ! Operate on value
298
    case('+')
299
      cval=oper1
300
    case('-')
301
      cval=-oper1
302
    end select
303
    call pushval(cval)             ! Put result back on stack
304
  case('F')  ! Token is a function name
305
    if(itop == 0) then
306
      ierr=7
307
      write(*,*) 'Error in expression ',trim(expr)
308
      write(*,*) 'Missing argument(s) for function ',trim(x%char)
309
      write(*,*)
310
      return
311
    else  
312
      call popval(oper1)           ! Pull top value off stack
313
    end if 
314
    tempstr=uppercase(x%char)
315
    select case(trim(tempstr))      ! Evaluate function
316
    case('SIN')
317
      cval=sin(oper1)
318
    case('COS')
319
      cval=cos(oper1)
320
    case('TAN')
321
      oper2=cos(oper1)
322
      if(abs(oper2) == 0.0_kr8) then
323
        ierr=14
324
        write(*,*) 'Error: argument of tan function a multiple',&
325
        ' of pi/2 in expression ',trim(expr)
326
        write(*,*)
327
        return
328
      else 
329
        cval=sin(oper1)/oper2
330
      endif
331
    case('SQRT')
332
      if(real(oper1,kr8) < 0. .and. aimag(oper1)==0.) then
333
        ierr=9
334
        write(*,*) 'Warning: square root of negative real number',&
335
                   ' in expression ',trim(expr)
336
        write(*,*)
337
      end if
338
      cval=sqrt(oper1)
339
    case('ABS')
340
      cval=abs(oper1)
341
    case('LN')
342
      if(real(oper1,kr8) <= 0. .and. aimag(oper1)==0.) then
343
        ierr=8
344
        write(*,*) 'Error: negative real or zero argument for',&
345
                   ' natural logarithm in expression ',trim(expr)
346
        write(*,*)
347
        return
348
      end if
349
      cval=log(oper1)
350
    case('LOG')
351
      if(real(oper1,kr8) <= 0. .and. aimag(oper1)==0.) then
352
        ierr=8
353
        write(*,*) 'Error: negative real or zero argument for base',&
354
                   '10 logarithm in expression ',trim(expr)
355
        write(*,*)
356
        return
357
      end if
358
      cval=log(oper1)/2.30258509299405_kr8
359
    case('EXP')
360
      cval=exp(oper1)
361
    case('COMPLEX')
362
      if(itop == 0) then
363
        ierr=7
364
        write(*,*) 'Error in expression ',trim(expr)
365
        write(*,*) 'Missing argument(s) for function ',trim(x%char)
366
        write(*,*)
367
        return
368
      else  
369
        call popval(oper2)  ! Pull second argument off stack
370
      end if 
371
      valr=real(oper2,kr8)
372
      vali=real(oper1,kr8)
373
      cval=cmplx(valr,vali,kc8)
374
    case('CONJG')
375
      cval=conjg(oper1)
376
    case('ANG')
377
      cval=atan2(aimag(oper1),real(oper1,kr8))
378
    case('REAL')
379
      cval=real(oper1,kr8)
380
    case('IMAG')
381
      cval=aimag(oper1)
382
    case default ! Undefined function
383
      ierr=13
384
      write(*,*) 'Error: the function ',trim(x%char), ' is undefined',&
385
                 ' in the expression ',trim(expr)
386
      write(*,*)
387
      return
388
    end select
389
    call pushval(cval)    ! Put result back on stack
390
  end select
391
end do
392
 
393
end subroutine evalexpr_dc
394
 
395
!**********************************************************************
396
 
397
SUBROUTINE GET_NEXT_TOKEN(str,tok,icp,isp)
398
 
399
character(len=*) :: str
400
character :: cop,chtemp
401
type(item) :: tok
402
integer :: icp
403
 
404
lstr=len_trim(str)
405
if(lstr == 0) then
406
  tok%char='#'             ! Output end token
407
  tok%type='E'
408
  return
409
end if
410
ipos=scan(str,'+-*/^(),')  ! Look for an arithmetic operator 
411
                           ! + - * / ^ ( ) or ,
412
cop=str(ipos:ipos)                 
413
select case (ipos)              
414
case(0)    ! Operators not present
415
  ntok=ntok+1
416
  tok%char=str
417
  tok%type='S'
418
  str=''
419
  icp=0
420
  isp=0
421
case(1) 
422
  tok%char=cop
423
  select case(cop)
424
  case('+','-')
425
    if(ibin==0) then
426
      tok%type='U'
427
      icp=4
428
      isp=3
429
    else
430
      tok%type='B'
431
      icp=1
432
      isp=1
433
    end if
434
    ibin=0
435
  case('*','/')
436
    tok%type='B'
437
    icp=2
438
    isp=2
439
    ibin=0
440
  case('^')
441
    tok%type='B'
442
    icp=4
443
    isp=3
444
    ibin=0
445
  case('(')
446
    tok%type='L'
447
    icp=4
448
    isp=0
449
    ibin=0
450
  case(')')
451
    tok%type='R'
452
    icp=0
453
    isp=0
454
    ibin=1
455
  case(',')
456
    tok%type='D'
457
    icp=0
458
    isp=0
459
    ibin=0
460
  end select
461
  str=str(2:)
462
case(2:)
463
  select case(cop)
464
  case('(')
465
    tok%char=str(1:ipos-1)
466
    tok%type='F'
467
    icp=4
468
    isp=0
469
    ibin=0
470
    str=str(ipos:)
471
  case('+','-')
472
    chtemp=uppercase(str(ipos-1:ipos-1))
473
    if(is_letter(str(1:1))==.true. .or. chtemp/='E') then
474
      tok%char=str(1:ipos-1)
475
      tok%type='S'
476
      icp=0
477
      isp=0
478
      ibin=1
479
      str=str(ipos:)
480
    else
481
      inext=scan(str(ipos+1:),'+-*/^(),')
482
      if(inext==0) then
483
        tok%char=str
484
        tok%type='S'
485
        icp=0
486
        isp=0
487
        ibin=0
488
        str=''
489
      else
490
        tok%char=str(1:ipos+inext-1)
491
        tok%type='S'
492
        icp=0
493
        isp=0
494
        ibin=1
495
        str=str(ipos+inext:)
496
      end if
497
    end if
498
  case default
499
    tok%char=str(1:ipos-1)
500
    tok%type='S'
501
    icp=0
502
    isp=0
503
    ibin=1
504
    str=str(ipos:)
505
  end select
506
end select
507
 
508
end subroutine get_next_token
509
 
510
 
511
!**********************************************************************
512
 
513
SUBROUTINE EVALEXPR_SC(expr,val)    ! Evaluate expression expr for
514
                                    ! val single precision complex
515
character(len=*) :: expr
516
complex(kc4) :: val
517
complex(kc8) :: vald
518
 
519
call evalexpr_dc(expr,vald)
520
val=vald
521
 
522
end subroutine evalexpr_sc
523
 
524
!**********************************************************************
525
 
526
SUBROUTINE EVALEXPR_SR(expr,val)    ! Evaluate expression expr for
527
                                    ! val single precision real
528
character(len=*) :: expr
529
real(kr4) :: val
530
complex(kc8) :: vald
531
 
532
call evalexpr_dc(expr,vald)
533
val=real(vald)
534
 
535
end subroutine evalexpr_sr
536
 
537
!**********************************************************************
538
 
539
SUBROUTINE EVALEXPR_DR(expr,val)    ! Evaluate expression expr for 
540
                                    ! val double precision real
541
character(len=*) :: expr
542
real(kr8) :: val
543
complex(kc8) :: vald
544
 
545
call evalexpr_dc(expr,vald)
546
val=real(vald,kr8)
547
 
548
end subroutine evalexpr_dr
549
 
550
!**********************************************************************
551
 
552
SUBROUTINE EVALEXPR_SI(expr,ival)   ! Evaluate expression expr for 
553
                                    ! ival single precision integer
554
character(len=*) :: expr
555
integer(ki4) :: ival
556
complex(kc8) :: vald
557
 
558
call evalexpr_dc(expr,vald)
559
ival=nint(real(vald,kr8),ki4)
560
 
561
end subroutine evalexpr_si
562
 
563
!**********************************************************************
564
 
565
SUBROUTINE EVALEXPR_DI(expr,ival)   ! Evaluate expression expr for 
566
                                    ! ival double precision integer
567
character(len=*) :: expr
568
integer(ki8) :: ival
569
complex(kc8) :: vald
570
 
571
call evalexpr_dc(expr,vald)
572
ival=nint(real(vald,kr8),ki8)
573
 
574
end subroutine evalexpr_di
575
 
576
 
577
!**********************************************************************
578
SUBROUTINE VALDEF_DC(sym,val)    ! Associates sym with val in symbol table,
579
                                 ! val double precision complex
580
character(len=*) :: sym
581
character(len=len_trim(sym)) :: usym
582
complex(kc8) :: val
583
 
584
ierr=0
585
if(nparams == 0) then               ! Initialize symbol table
586
  params(1)%symbol='PI'
587
  params(1)%value=(3.14159265358979_kr8,0.0_kr8)
588
  params(2)%symbol='I'
589
  params(2)%value=(0.0_kr8,1.0_kr8)
590
  nparams=2
591
end if
592
 
593
! Assign val to sym if sym is already in symbol table
594
usym=uppercase(sym)
595
if(is_letter(sym(1:1))==.false. .or. len_trim(sym)>24) then
596
  ierr=11
597
  write(*,*) 'Error: symbol ',trim(sym),' has improper format'
598
  write(*,*)
599
  return
600
end if
601
do i=1,nparams
602
  if(trim(usym)==trim(params(i)%symbol)) then
603
    params(i)%value=val
604
    return
605
  end if
606
end do
607
 
608
nparams=nparams+1    ! Otherwise assign val to new symbol sym
609
params(nparams)%symbol=usym
610
params(nparams)%value=val
611
 
612
end subroutine valdef_dc
613
 
614
 
615
!**********************************************************************
616
 
617
SUBROUTINE VALDEF_SC(sym,val)     ! Associates sym with val in symbol table,
618
                                  ! val single precision complex
619
character(len=*) :: sym
620
complex(kc4) :: val
621
complex(kc8) :: vald
622
 
623
vald=val
624
call valdef_dc(sym,vald)
625
 
626
end subroutine valdef_sc
627
 
628
 
629
!**********************************************************************
630
 
631
SUBROUTINE VALDEF_DR(sym,val)    ! Associates sym with val in symbol table,
632
                                 ! val double precision real
633
character(len=*) :: sym
634
real(kr8) :: val
635
complex(kc8) :: vald
636
 
637
vald=cmplx(val,0.0_kr8,kc8)
638
call valdef_dc(sym,vald)
639
 
640
end subroutine valdef_dr
641
 
642
 
643
!**********************************************************************
644
 
645
SUBROUTINE VALDEF_SR(sym,val)    ! Associates sym with val in symbol table,
646
                                 ! val single precision real
647
character(len=*) :: sym
648
real(kr4) :: val
649
complex(kc8) :: vald
650
 
651
vald=cmplx(val,0.0,kc8)
652
call valdef_dc(sym,vald)
653
 
654
end subroutine valdef_sr
655
 
656
 
657
!**********************************************************************
658
 
659
SUBROUTINE VALDEF_DI(sym,ival)   ! Associates sym with ival in symbol table,
660
                                 ! ival double precision integer 
661
character(len=*) :: sym
662
integer(ki8) :: ival
663
complex(kc8) :: vald
664
 
665
vald=cmplx(real(ival,kr8),0.0_kr8,kc8)
666
call valdef_dc(sym,vald)
667
 
668
end subroutine valdef_di
669
 
670
 
671
!**********************************************************************
672
 
673
SUBROUTINE VALDEF_SI(sym,ival)   ! Associates sym with ival in symbol table,
674
                                 ! ival single precision integer
675
character(len=*) :: sym
676
integer(ki4) :: ival
677
complex(kc8) :: vald
678
 
679
vald=cmplx(real(ival,kr8),0.0,kc8)
680
call valdef_dc(sym,vald)
681
 
682
end subroutine valdef_si
683
 
684
 
685
!**********************************************************************
686
 
687
SUBROUTINE STRDEF(sym,expr)      ! Associates sym with the value of the
688
                                 ! expression expr
689
 
690
character(len=*) :: sym,expr
691
complex(kc8) :: val
692
 
693
if(nparams == 0) then            ! Initialize symbol table
694
  params(1)%symbol='PI'
695
  params(1)%value=(3.14159265358979_kr8,0.0_kr8)
696
  params(2)%symbol='I'
697
  params(2)%value=(0.0_kr8,1.0_kr8)
698
  nparams=2
699
end if
700
 
701
call evalexpr_dc(expr,val)       ! val is value of expression expr
702
if(ierr==0 .or. ierr==9) then
703
  call valdef_dc(sym,val)          ! Assign val to symbol sym
704
end if
705
 
706
end subroutine strdef
707
 
708
 
709
!**********************************************************************
710
 
711
SUBROUTINE VALUEP(xinchar,cval)  ! Finds double precision complex value 
712
                                 ! corresponding to number string xinchar 
713
                                 ! or value in symbol table corresponding 
714
                                 ! to symbol name xinchar.
715
 
716
character (len=*):: xinchar
717
complex(kc8) :: cval
718
real(kr8) :: rval
719
 
720
ierr=0
721
 
722
if(is_letter(xinchar(1:1))==.true.) then   ! xinchar is a symbol
723
  call getparam(xinchar,cval)
724
else                               ! xinchar is a number string
725
  call value(xinchar,rval,ios)     ! rval is the value of xinchar
726
  if(ios > 0) then
727
    ierr=3
728
    write(*,*) 'Error: number string ',trim(xinchar),' does not correspond to a valid number' 
729
    write(*,*)
730
  end if
731
  cval=cmplx(rval,0.0_kr8,kc8)
732
  return
733
end if
734
 
735
end subroutine valuep
736
 
737
 
738
!**********************************************************************
739
 
740
 
741
SUBROUTINE PUSHOP(op)  ! Puts an operator on operator stack
742
 
743
type(item):: op
744
 
745
itop=itop+1
746
if(itop > numtok) then
747
  write(*,*) 'Error: operator stack overflow in evaluation of expression'
748
  write(*,*)
749
  return
750
end if
751
opstack(itop)=op
752
 
753
end subroutine pushop
754
 
755
SUBROUTINE POPOP(op) ! Takes top operator of operator stack and assigns it to op
756
 
757
type(item):: op
758
 
759
op=opstack(itop)
760
itop=itop-1
761
 
762
end subroutine popop
763
 
764
SUBROUTINE PUSHVAL(val) ! Puts value on value stack
765
 
766
complex(kc8) :: val
767
 
768
itop=itop+1
769
if(itop > numtok) then
770
  write(*,*) 'Error: value stack overflow in evaluation of expression'
771
  write(*,*)
772
  return
773
end if
774
valstack(itop)=val
775
 
776
end subroutine pushval
777
 
778
SUBROUTINE POPVAL(val) ! Takes top value off value stack and assigns it to val
779
 
780
complex(kc8) :: val
781
 
782
val=valstack(itop)
783
itop=itop-1
784
 
785
end subroutine popval
786
 
787
!**********************************************************************
788
 
789
SUBROUTINE GETPARAM_DC(sym,var)  ! Find double precision complex value var
790
                                 ! corresponding to symbol sym
791
 
792
character(len=*) :: sym
793
character(len=len_trim(sym)) :: usym
794
complex(kc8) :: var
795
 
796
ierr=0
797
sym=adjustl(sym)
798
if(is_letter(sym(1:1))==.false. .or. len_trim(sym)>24) then
799
  ierr=11
800
  write(*,*) 'Error: symbol ',trim(sym),' has incorrect format'
801
  write(*,*)
802
  return
803
end if
804
ifind=0
805
usym=uppercase(sym)
806
do j=1,nparams
807
  if(trim(usym) == trim(params(j)%symbol)) then
808
    var=params(j)%value
809
    ifind=j
810
    exit
811
  end if
812
end do
813
if(ifind == 0) then          
814
  ierr=4
815
  write(*,*) 'Error: symbol ',trim(sym), ' not in symbol table'
816
  write(*,*) 
817
  return
818
end if
819
 
820
end subroutine getparam_dc
821
 
822
!**********************************************************************
823
 
824
SUBROUTINE GETPARAM_SC(sym,var)  ! Find single precision complex value var
825
                                 ! corresponding to symbol sym
826
 
827
 
828
character(len=*) :: sym
829
complex(kc4) :: var
830
complex(kc8) :: vard
831
 
832
call getparam_dc(sym,vard)
833
var=vard
834
 
835
end subroutine getparam_sc
836
 
837
!**********************************************************************
838
 
839
SUBROUTINE GETPARAM_DR(sym,var)  ! Find double precision real value var
840
                                 ! corresponding to symbol sym
841
 
842
 
843
character(len=*) :: sym
844
real(kr8) :: var
845
complex(kc8) :: vard
846
 
847
call getparam_dc(sym,vard)
848
var=real(vard,kr8)
849
 
850
end subroutine getparam_dr
851
 
852
!**********************************************************************
853
 
854
SUBROUTINE GETPARAM_SR(sym,var)  ! Find single precision real value var
855
                                 ! corresponding to symbol sym
856
 
857
 
858
character(len=*) :: sym
859
real(kr4) :: var
860
complex(kc8) :: vard
861
 
862
call getparam_dc(sym,vard)
863
var=real(vard)
864
 
865
end subroutine getparam_sr
866
 
867
!**********************************************************************
868
 
869
SUBROUTINE GETPARAM_DI(sym,ivar)  ! Find double precision integer value ivar
870
                                  ! corresponding to symbol sym
871
 
872
 
873
character(len=*) :: sym
874
integer(ki8) :: ivar
875
complex(kc8) :: vard
876
 
877
call getparam_dc(sym,vard)
878
ivar=nint(real(vard,kr8),ki8)
879
 
880
end subroutine getparam_di
881
 
882
!**********************************************************************
883
 
884
SUBROUTINE GETPARAM_SI(sym,ivar)  ! Find single precision integer value ivar
885
                                  ! corresponding to symbol sym
886
 
887
 
888
character(len=*) :: sym
889
integer(ki4) :: ivar
890
complex(kc8) :: vard
891
 
892
call getparam_dc(sym,vard)
893
ivar=nint(real(vard,kr8),ki4)
894
 
895
end subroutine getparam_si
896
 
897
!**********************************************************************
898
 
899
SUBROUTINE EVALEQN(eqn)  ! Evaluate an equation
900
 
901
character(len=*) :: eqn
902
character(len=len(eqn)) :: args(2)
903
complex(kc8) :: val
904
 
905
call parse(eqn,'=',args,nargs)   ! Seperate right- and left-hand-sides
906
call defparam(adjustl(args(1)),args(2)) ! Evaluate right-hand-side and
907
                                        ! assign to symbol on the
908
                                        ! left-hand-side.
909
end subroutine evaleqn
910
 
911
!**********************************************************************
912
 
913
SUBROUTINE LISTVAR      ! List all variables and their values
914
 
915
write(*,'(/a)') ' VARIABLE LIST:'
916
if(nparams == 0) then            ! Initialize symbol table
917
  params(1)%symbol='PI'
918
  params(1)%value=(3.14159265358979_kr8,0.0_kr8)
919
  params(2)%symbol='I'
920
  params(2)%value=(0.0_kr8,1.0_kr8)
921
  nparams=2
922
end if
923
do i=1,nparams
924
  write(*,*) trim(params(i)%symbol),' = ',params(i)%value
925
end do
926
 
927
end subroutine listvar
928
 
929
!**********************************************************************
930
 
931
end module evaluate