trunk/IOIPSL/mathelp.f

MODULE mathelp

  ! From mathelp.f90, version 2.0 2004/04/05 14:47:50

  implicit none

  PRIVATE
  PUBLIC :: moycum, trans_buff, buildop

  !- Variables used to detect and identify the operations
  CHARACTER(LEN=80), SAVE :: seps='( ) , + - / * ^'

CONTAINS

  SUBROUTINE buildop (str, ex_topps, topp, nbops_max, &
       &                  missing_val, opps, scal, nbops)
    !- This subroutine decomposes the input string in the elementary
    !- functions which need to be applied to the vector of data.
    !- This vector is represented by X in the string.
    !- This subroutine is the driver of the decomposition and gets
    !- the time operation but then call decoop for the other operations
    !- INPUT

    !- str      : String containing the operations
    !- ex_toops : The time operations that can be expected
    !-            within the string

    !- OUTPUT

    USE errioipsl, ONLY : histerr

    CHARACTER(LEN=80) :: str
    CHARACTER(LEN=*) :: ex_topps
    CHARACTER(LEN=7) :: topp
    INTEGER :: nbops_max, nbops
    CHARACTER(LEN=7) :: opps(nbops_max)
    REAL :: scal(nbops_max), missing_val

    CHARACTER(LEN=80) :: new_str
    INTEGER :: leng, ind_opb, ind_clb

    LOGICAL :: check = .FALSE.
    !---------------------------------------------------------------------
    IF (check) WRITE(*, *) 'buildop : Some preliminary cleaning'

    leng = LEN_TRIM(str)
    IF ( str(1:1) == '(' .AND. str(leng:leng) == ')' ) THEN
       str = str(2:leng-1)
       leng = leng-2
    ENDIF

    IF (check) &
         &  WRITE(*, *) 'buildop : Starting to test the various options'

    IF (leng <= 5 .AND. INDEX(ex_topps, str(1:leng)) > 0) THEN
       IF (check) WRITE(*, *) 'buildop : Time operation only'
       nbops = 0
       topp = str(1:leng)
    ELSE
       IF (check) THEN
          WRITE(*, *) 'buildop : Time operation and something else'
       ENDIF
       !--
       ind_opb = INDEX(str(1:leng), '(')
       IF (ind_opb > 0) THEN
          IF (INDEX(ex_topps, str(1:ind_opb-1)) > 0) THEN
             IF (check) THEN
                WRITE(*, '(2a)') &
                     &          ' buildop : Extract time operation from : ', str
             ENDIF
             topp = str(1:ind_opb-1)
             ind_clb = INDEX(str(1:leng), ')', BACK=.TRUE.)
             new_str = str(ind_opb+1:ind_clb-1)
             IF (check) THEN
                WRITE(*, '(2a, 2I3)') &
                     &          ' buildop : Call decoop ', new_str, ind_opb, ind_clb
             ENDIF
             CALL decoop (new_str, nbops_max, missing_val, opps, scal, nbops)
          ELSE
             CALL histerr(3, 'buildop', &
                  &        'time opperation does not exist', str(1:ind_opb-1), ' ')
          ENDIF
       ELSE
          CALL histerr(3, 'buildop', &
               &      'some long opperation exists but wihout parenthesis', &
               &      str(1:leng), ' ')
       ENDIF
    ENDIF

    IF (check) THEN
       DO leng=1, nbops
          WRITE(*, *) &
               &      'buildop : i -- opps, scal : ', leng, opps(leng), scal(leng)
       ENDDO
    ENDIF
    !---------------------
  END SUBROUTINE buildop

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

  SUBROUTINE decoop (pstr, nbops_max, missing_val, opps, scal, nbops)
    USE errioipsl, ONLY : histerr

    CHARACTER(LEN=80) :: pstr
    INTEGER :: nbops_max, nbops
    CHARACTER(LEN=7) :: opps(nbops_max)
    REAL :: scal(nbops_max), missing_val

    CHARACTER(LEN=1) :: f_char(2), s_char(2)
    INTEGER :: nbsep, f_pos(2), s_pos(2)
    CHARACTER(LEN=20) :: opp_str, scal_str
    CHARACTER(LEN=80) :: str
    INTEGER :: xpos, leng, ppos, epos, int_tmp
    CHARACTER(LEN=3) :: tl, dl
    CHARACTER(LEN=10) :: fmt

    LOGICAL :: check = .FALSE., prio
    CHARACTER(LEN=80), SAVE :: ops = '+ - * / ^'
    CHARACTER(LEN=80), SAVE :: mima = 'min max'
    CHARACTER(LEN=250), SAVE :: funcs = &
         'sin cos tan asin acos atan exp log sqrt chs abs ' &
         //'cels kelv deg rad gather scatter fill coll undef only ident'
    !---------------------------------------------------------------------
    IF (check) WRITE(*, '(2a)') ' decoop : Incoming string : ', pstr

    nbops = 0
    str = pstr

    CALL findsep (str, nbsep, f_char, f_pos, s_char, s_pos)
    IF (check) WRITE(*, *) 'decoop : Out of findsep', nbsep
    DO WHILE (nbsep > 0)
       xpos = INDEX(str, 'X')
       leng = LEN_TRIM(str)
       nbops = nbops+1
       !--
       IF (check) THEN
          WRITE(*, *) 'decoop : str -->', str(1:leng)
          WRITE(*, *) s_char(1), '-', f_char(1), '|', f_char(2), '-', s_char(2)
          WRITE(*, *) s_pos(1), '-', f_pos(1), '|', f_pos(2), '-', s_pos(2)
       ENDIF
       !--
       IF (nbops > nbops_max-1) THEN
          CALL histerr(3, 'decoop', 'Expression too complex', str, ' ')
       ENDIF
       !--
       IF (check) WRITE(*, *) 'decoop : --', nbops, ' ', str(1:leng)
       !---
       !-- Start the analysis of the syntax. 3 types of constructs
       !-- are recognized.  They are scanned sequentialy
       !---
       IF (nbsep == 1) THEN
          IF (check) WRITE(*, *) 'decoop : Only one operation'
          IF (INDEX(ops, f_char(1)) > 0) THEN
             !------ Type : scal+X
             IF (f_char(1) == '-' .OR. f_char(1) == '/') THEN
                opp_str = f_char(1)//'I'
             ELSE
                opp_str = f_char(1)
             ENDIF
             scal_str = str(s_pos(1)+1:f_pos(1)-1)
             str = 'X'
          ELSE IF (INDEX(ops, f_char(2)) > 0) THEN
             !------ Type : X+scal
             opp_str = f_char(2)
             scal_str = str(f_pos(2)+1:s_pos(2)-1)
             str = 'X'
          ELSE
             CALL histerr(3, 'decoop', &
                  &        'Unknown operations of type X+scal', f_char(1), pstr)
          ENDIF
       ELSE
          IF (check) WRITE(*, *) 'decoop : More complex operation'
          IF ( f_char(1) == '(' .AND. f_char(2) == ')' ) THEN
             !------ Type : sin(X)
             opp_str = str(s_pos(1)+1:f_pos(1)-1)
             scal_str = '?'
             str = str(1:s_pos(1))//'X'//str(f_pos(2)+1:leng)
          ELSE IF (    (f_char(1) == '(' .AND. f_char(2) == ', ')&
               &             .OR.(f_char(1) == ', ' .AND. f_char(2) == ')')) THEN
             !------ Type : max(X, scal) or max(scal, X)
             IF (f_char(1) == '(' .AND. s_char(2) == ')') THEN
                !-------- Type : max(X, scal)
                opp_str = str(f_pos(1)-3:f_pos(1)-1)
                scal_str = str(f_pos(2)+1:s_pos(2)-1)
                str = str(1:f_pos(1)-4)//'X'//str(s_pos(2)+1:leng)
             ELSE IF (f_char(1) == ', ' .AND. s_char(1) == '(') THEN
                !-------- Type : max(scal, X)
                opp_str = str(s_pos(1)-3:s_pos(1)-1)
                scal_str = str(s_pos(1)+1:f_pos(1)-1)
                str = str(1:s_pos(1)-4)//'X'//str(f_pos(2)+1:leng)
             ELSE
                CALL histerr(3, 'decoop', 'Syntax error 1', str, ' ')
             ENDIF
          ELSE
             prio = (f_char(2) == '*').OR.(f_char(2) == '^')
             IF (     (INDEX(ops, f_char(1)) > 0) &
                  &          .AND.(xpos-f_pos(1) == 1).AND.(.NOT.prio) ) THEN
                !-------- Type : ... scal+X ...
                IF (f_char(1) == '-' .OR. f_char(1) == '/') THEN
                   opp_str = f_char(1)//'I'
                ELSE
                   opp_str = f_char(1)
                ENDIF
                scal_str = str(s_pos(1)+1:f_pos(1)-1)
                str = str(1:s_pos(1))//'X'//str(f_pos(1)+2:leng)
             ELSE IF (     (INDEX(ops, f_char(2)) > 0) &
                  &               .AND.(f_pos(2)-xpos == 1) ) THEN
                !-------- Type : ... X+scal ...
                opp_str = f_char(2)
                scal_str = str(f_pos(2)+1:s_pos(2)-1)
                str = str(1:f_pos(2)-2)//'X'//str(s_pos(2):leng)
             ELSE
                CALL histerr(3, 'decoop', 'Syntax error 2', str, ' ')
             ENDIF
          ENDIF
       ENDIF
       !---
       IF (check) WRITE(*, *) 'decoop : Finished syntax, str = ', TRIM(str)
       !---
       !-- Now that the different components of the operation are identified
       !-- we transform them into what is going to be used in the program
       !---
       IF (INDEX(scal_str, '?') > 0) THEN
          IF (INDEX(funcs, opp_str(1:LEN_TRIM(opp_str))) > 0) THEN
             opps(nbops) = opp_str(1:LEN_TRIM(opp_str))
             scal(nbops) =  missing_val
          ELSE
             CALL histerr(3, 'decoop', &
                  &        'Unknown function', opp_str(1:LEN_TRIM(opp_str)), ' ')
          ENDIF
       ELSE
          leng = LEN_TRIM(opp_str)
          IF (INDEX(mima, opp_str(1:leng)) > 0) THEN
             opps(nbops) = 'fu'//opp_str(1:leng)
          ELSE
             IF (INDEX(opp_str(1:leng), '+') > 0) THEN
                opps(nbops) = 'add'
             ELSE IF (INDEX(opp_str(1:leng), '-I') > 0) THEN
                opps(nbops) = 'subi'
             ELSE IF (INDEX(opp_str(1:leng), '-') > 0) THEN
                opps(nbops) = 'sub'
             ELSE IF (INDEX(opp_str(1:leng), '*') > 0) THEN
                opps(nbops) = 'mult'
             ELSE IF (INDEX(opp_str(1:leng), '/') > 0) THEN
                opps(nbops) = 'div'
             ELSE IF (INDEX(opp_str(1:leng), '/I') > 0) THEN
                opps(nbops) = 'divi'
             ELSE IF (INDEX(opp_str(1:leng), '^') > 0) THEN
                opps(nbops) = 'power'
             ELSE
                CALL histerr(3, 'decoop', &
                     &          'Unknown operation', opp_str(1:leng), ' ')
             ENDIF
          ENDIF
          !-----
          leng = LEN_TRIM(scal_str)
          ppos = INDEX(scal_str, '.')
          epos = INDEX(scal_str, 'e')
          IF (epos == 0) epos = INDEX(scal_str, 'E')
          !-----
          !---- Try to catch a few errors
          !-----
          IF (INDEX(ops, scal_str) > 0) THEN
             CALL histerr(3, 'decoop', &
                  &        'Strange scalar you have here ', scal_str, pstr)
          ENDIF
          IF (epos > 0) THEN
             WRITE(tl, '(I3.3)') leng
             WRITE(dl, '(I3.3)') epos-ppos-1
             fmt='(e'//tl//'.'//dl//')'
             READ(scal_str, fmt) scal(nbops)
          ELSE IF (ppos > 0) THEN
             WRITE(tl, '(I3.3)') leng
             WRITE(dl, '(I3.3)') leng-ppos
             fmt='(f'//tl//'.'//dl//')'
             READ(scal_str, fmt) scal(nbops)
          ELSE
             WRITE(tl, '(I3.3)') leng
             fmt = '(I'//tl//')'
             READ(scal_str, fmt) int_tmp
             scal(nbops) = REAL(int_tmp)
          ENDIF
       ENDIF
       IF (check) WRITE(*, *) 'decoop : Finished interpretation'
       CALL findsep(str, nbsep, f_char, f_pos, s_char, s_pos)
    ENDDO
    !--------------------
  END SUBROUTINE decoop

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

  SUBROUTINE findsep (str, nbsep, f_char, f_pos, s_char, s_pos)
    !- Subroutine finds all separators in a given string
    !- It returns the following information about str :
    !-   f_char : The first separation character
    !-            (1 for before and 2 for after)
    !-   f_pos  : The position of the first separator
    !-   s_char : The second separation character
    !-            (1 for before and 2 for after)
    !-   s_pos  : The position of the second separator
    USE errioipsl, ONLY : histerr

    CHARACTER(LEN=80) :: str
    INTEGER :: nbsep
    CHARACTER(LEN=1), DIMENSION(2) :: f_char, s_char
    INTEGER, DIMENSION(2) :: f_pos, s_pos

    CHARACTER(LEN=70) :: str_tmp
    LOGICAL :: f_found, s_found
    INTEGER :: ind, xpos, leng, i

    LOGICAL :: check = .FALSE.
    !---------------------------------------------------------------------
    IF (check) WRITE(*, *) 'findsep : call cleanstr: ', TRIM(str)

    CALL cleanstr(str)

    IF (check) WRITE(*, *) 'findsep : out of cleanstr: ', TRIM(str)

    xpos = INDEX(str, 'X')
    leng = LEN_TRIM(str)

    f_pos(1:2) = (/ 0, leng+1 /)
    f_char(1:2) = (/ '?', '?' /)
    s_pos(1:2) = (/ 0, leng+1 /)
    s_char(1:2) = (/ '?', '?' /)

    nbsep = 0

    f_found = .FALSE.
    s_found = .FALSE.
    IF (xpos > 1) THEN
       DO i=xpos-1, 1, -1
          ind = INDEX(seps, str(i:i))
          IF (ind > 0) THEN
             IF (.NOT.f_found) THEN
                f_char(1) = str(i:i)
                f_pos(1) = i
                nbsep = nbsep+1
                f_found = .TRUE.
             ELSE IF (.NOT.s_found) THEN
                s_char(1) = str(i:i)
                s_pos(1) = i
                nbsep = nbsep+1
                s_found = .TRUE.
             ENDIF
          ENDIF
       ENDDO
    ENDIF

    f_found = .FALSE.
    s_found = .FALSE.
    IF (xpos < leng) THEN
       DO i=xpos+1, leng
          ind = INDEX(seps, str(i:i))
          IF (ind > 0) THEN
             IF (.NOT.f_found) THEN
                f_char(2) = str(i:i)
                f_pos(2) = i
                nbsep = nbsep+1
                f_found = .TRUE.
             ELSE IF (.NOT.s_found) THEN
                s_char(2) = str(i:i)
                s_pos(2) = i
                nbsep = nbsep+1
                s_found = .TRUE.
             ENDIF
          ENDIF
       ENDDO
    ENDIF

    IF (nbsep > 4) THEN
       WRITE(str_tmp, '("number :", I3)') nbsep
       CALL histerr(3, 'findsep', &
            &    'How can I find that many separators', str_tmp, str)
    ENDIF

    IF (check) WRITE(*, *) 'Finished findsep : ', nbsep, leng
    !---------------------
  END SUBROUTINE findsep

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

  SUBROUTINE cleanstr(str)
    !- We clean up the string by taking out the extra () and puting
    !- everything in lower case except for the X describing the variable
    use strlowercase_m, only: strlowercase

    CHARACTER(LEN=80) :: str

    INTEGER :: ind, leng, ic, it
    LOGICAL :: check = .FALSE.
    !---------------------------------------------------------------------
    leng = LEN_TRIM(str)
    CALL strlowercase(str)

    ind = INDEX(str, 'x')
    IF (check) THEN
       WRITE (*, *) 'cleanstr 1.0 : ind = ', ind, &
            &               ' str = ', str(1:leng), '---'
    ENDIF

    ! If the character before the x is not a letter then we can assume
    ! that it is the variable and promote it to a capital letter

    DO WHILE (ind > 0)
       ic = 0
       IF (ind > 1) ic = IACHAR(str(ind-1:ind-1))
       IF (ic < 97 .OR. ic > 122) THEN
          str(ind:ind) = 'X'
       ENDIF
       it = INDEX(str(ind+1:leng), 'x')
       IF (it > 0) THEN
          ind = ind+it
       ELSE
          ind = it
       ENDIF
    ENDDO

    IF (check) WRITE (*, *) 'cleanstr 2.0 : str = ', str(1:leng), '---'

    IF ( str(1:1) == '(' .AND. str(leng:leng) == ')' ) THEN
       str = str(2:leng-1)
    ENDIF

    IF (check) WRITE (*, *) 'cleanstr 3.0 : str = ', str(1:leng), '---'

    leng = LEN_TRIM(str)
    ind = INDEX(str, '((X))')
    IF (ind > 0) THEN
       str=str(1:ind-1)//'(X)'//str(ind+5:leng)//'  '
    ENDIF

    IF (check) WRITE (*, *) 'cleanstr 4.0 : str = ', str(1:leng), '---'

    leng = LEN_TRIM(str)
    ind = INDEX(str, '(X)')
    IF (ind > 0 .AND. ind+3 < leng) THEN
       IF (      (INDEX(seps, str(ind-1:ind-1)) > 0) &
            &      .AND. (INDEX(seps, str(ind+3:ind+3)) > 0) ) THEN
          str=str(1:ind-1)//'X'//str(ind+3:leng)//'  '
       ENDIF
    ENDIF

    IF (check) WRITE (*, *) 'cleanstr 5.0 : str = ', str(1:leng), '---'

    leng = LEN_TRIM(str)
    ind = INDEX(str(1:leng), ' ')
    DO WHILE (ind > 0)
       str=str(1:ind-1)//str(ind+1:leng)//' '
       leng = LEN_TRIM(str)
       ind = INDEX(str(1:leng), ' ')
    ENDDO

    IF (check) WRITE (*, *) 'cleanstr 6.0 : str = ', str(1:leng), '---'
    !----------------------
  END SUBROUTINE cleanstr

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

  SUBROUTINE moycum (opp, np, px, py, pwx)
    !- Does time operations
    USE errioipsl, ONLY : histerr

    CHARACTER(LEN=7) :: opp
    INTEGER :: np
    REAL, DIMENSION(:) :: px, py
    INTEGER :: pwx
    !---------------------------------------------------------------------
    IF (pwx /= 0) THEN
       IF      (opp == 'ave') THEN
          px(1:np)=(px(1:np)*pwx+py(1:np))/REAL(pwx+1)
       ELSE IF (opp == 't_sum') THEN
          px(1:np)=px(1:np)+py(1:np)
       ELSE IF ( (opp == 'l_min').OR.(opp == 't_min') ) THEN
          px(1:np)=MIN(px(1:np), py(1:np))
       ELSE IF ( (opp == 'l_max').OR.(opp == 't_max') ) THEN
          px(1:np)=MAX(px(1:np), py(1:np))
       ELSE
          CALL histerr(3, "moycum", 'Unknown time operation', opp, ' ')
       ENDIF
    ELSE
       IF      (opp == 'l_min') THEN
          px(1:np)=MIN(px(1:np), py(1:np))
       ELSE IF (opp == 'l_max') THEN
          px(1:np)=MAX(px(1:np), py(1:np))
       ELSE
          px(1:np)=py(1:np)
       ENDIF
    ENDIF
    !--------------------
  END SUBROUTINE moycum

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

  SUBROUTINE trans_buff (ox, sx, oy, sy, oz, sz, xsz, ysz, zsz, v3d, sl, v1d)
    !- This subroutine extracts from the full 3D variable the slab of
    !- data that will be used later. Perhaps there are hardware routines
    !- for this task on some computers. This routine will be obsolete in
    !- a F90 environnement

    !- INPUT
    !- ox  : Origin of slab of data in X
    !- sx  : Size of slab in X
    !- oy  : Origin of slab of data in Y
    !- sy  : Size of slab in Y
    !- oz  : Origin of slab of data in Z
    !- sz  : Size of slab in Z
    !- xsz, ysz, zsz : 3 sizes of full variable v3d
    !- v3d : The full 3D variable
    !- sl  : size of variable v1d
    !- v1d : The 1D variable containing the slab

    INTEGER :: ox, sx, oy, sy, oz, sz
    INTEGER :: xsz, ysz, zsz
    INTEGER :: sl
    REAL :: v3d(xsz, ysz, zsz)
    REAL :: v1d(sl)

    !---------------------------------------------------------------------

    V1d(:sx*sy*sz) = reshape(V3d(ox:ox-1+sx, oy:oy-1+sy, oz:oz-1+sz), &
         SHAPE=(/sx*sy*sz/))

  END SUBROUTINE trans_buff

END MODULE mathelp
1	guez	30	MODULE mathelp
2
3			! From mathelp.f90, version 2.0 2004/04/05 14:47:50
4
5			implicit none
6
7			PRIVATE
8	guez	32	PUBLIC :: moycum, trans_buff, buildop
9	guez	30
10			!- Variables used to detect and identify the operations
11	guez	32	CHARACTER(LEN=80), SAVE :: seps='( ) , + - / * ^'
12	guez	30
13			CONTAINS
14
15			SUBROUTINE buildop (str, ex_topps, topp, nbops_max, &
16			& missing_val, opps, scal, nbops)
17			!- This subroutine decomposes the input string in the elementary
18			!- functions which need to be applied to the vector of data.
19			!- This vector is represented by X in the string.
20			!- This subroutine is the driver of the decomposition and gets
21			!- the time operation but then call decoop for the other operations
22			!- INPUT
23
24			!- str : String containing the operations
25			!- ex_toops : The time operations that can be expected
26			!- within the string
27
28			!- OUTPUT
29
30	guez	32	USE errioipsl, ONLY : histerr
31
32	guez	30	CHARACTER(LEN=80) :: str
33			CHARACTER(LEN=*) :: ex_topps
34			CHARACTER(LEN=7) :: topp
35			INTEGER :: nbops_max, nbops
36			CHARACTER(LEN=7) :: opps(nbops_max)
37			REAL :: scal(nbops_max), missing_val
38
39			CHARACTER(LEN=80) :: new_str
40			INTEGER :: leng, ind_opb, ind_clb
41
42			LOGICAL :: check = .FALSE.
43			!---------------------------------------------------------------------
44			IF (check) WRITE(, ) 'buildop : Some preliminary cleaning'
45
46			leng = LEN_TRIM(str)
47			IF ( str(1:1) == '(' .AND. str(leng:leng) == ')' ) THEN
48			str = str(2:leng-1)
49			leng = leng-2
50			ENDIF
51
52			IF (check) &
53			& WRITE(, ) 'buildop : Starting to test the various options'
54
55			IF (leng <= 5 .AND. INDEX(ex_topps, str(1:leng)) > 0) THEN
56			IF (check) WRITE(, ) 'buildop : Time operation only'
57			nbops = 0
58			topp = str(1:leng)
59			ELSE
60			IF (check) THEN
61			WRITE(, ) 'buildop : Time operation and something else'
62			ENDIF
63			!--
64			ind_opb = INDEX(str(1:leng), '(')
65			IF (ind_opb > 0) THEN
66			IF (INDEX(ex_topps, str(1:ind_opb-1)) > 0) THEN
67			IF (check) THEN
68			WRITE(*, '(2a)') &
69			& ' buildop : Extract time operation from : ', str
70			ENDIF
71			topp = str(1:ind_opb-1)
72			ind_clb = INDEX(str(1:leng), ')', BACK=.TRUE.)
73			new_str = str(ind_opb+1:ind_clb-1)
74			IF (check) THEN
75			WRITE(*, '(2a, 2I3)') &
76			& ' buildop : Call decoop ', new_str, ind_opb, ind_clb
77			ENDIF
78			CALL decoop (new_str, nbops_max, missing_val, opps, scal, nbops)
79			ELSE
80			CALL histerr(3, 'buildop', &
81			& 'time opperation does not exist', str(1:ind_opb-1), ' ')
82			ENDIF
83			ELSE
84			CALL histerr(3, 'buildop', &
85			& 'some long opperation exists but wihout parenthesis', &
86			& str(1:leng), ' ')
87			ENDIF
88			ENDIF
89
90			IF (check) THEN
91			DO leng=1, nbops
92			WRITE(, ) &
93			& 'buildop : i -- opps, scal : ', leng, opps(leng), scal(leng)
94			ENDDO
95			ENDIF
96			!---------------------
97			END SUBROUTINE buildop
98
99			!************************************************
100
101			SUBROUTINE decoop (pstr, nbops_max, missing_val, opps, scal, nbops)
102	guez	32	USE errioipsl, ONLY : histerr
103
104	guez	30	CHARACTER(LEN=80) :: pstr
105			INTEGER :: nbops_max, nbops
106			CHARACTER(LEN=7) :: opps(nbops_max)
107			REAL :: scal(nbops_max), missing_val
108
109			CHARACTER(LEN=1) :: f_char(2), s_char(2)
110			INTEGER :: nbsep, f_pos(2), s_pos(2)
111			CHARACTER(LEN=20) :: opp_str, scal_str
112			CHARACTER(LEN=80) :: str
113			INTEGER :: xpos, leng, ppos, epos, int_tmp
114			CHARACTER(LEN=3) :: tl, dl
115			CHARACTER(LEN=10) :: fmt
116
117			LOGICAL :: check = .FALSE., prio
118	guez	32	CHARACTER(LEN=80), SAVE :: ops = '+ - * / ^'
119			CHARACTER(LEN=80), SAVE :: mima = 'min max'
120			CHARACTER(LEN=250), SAVE :: funcs = &
121			'sin cos tan asin acos atan exp log sqrt chs abs ' &
122			//'cels kelv deg rad gather scatter fill coll undef only ident'
123	guez	30	!---------------------------------------------------------------------
124			IF (check) WRITE(*, '(2a)') ' decoop : Incoming string : ', pstr
125
126			nbops = 0
127			str = pstr
128
129			CALL findsep (str, nbsep, f_char, f_pos, s_char, s_pos)
130			IF (check) WRITE(, ) 'decoop : Out of findsep', nbsep
131			DO WHILE (nbsep > 0)
132			xpos = INDEX(str, 'X')
133			leng = LEN_TRIM(str)
134			nbops = nbops+1
135			!--
136			IF (check) THEN
137			WRITE(, ) 'decoop : str -->', str(1:leng)
138			WRITE(, ) s_char(1), '-', f_char(1), '\|', f_char(2), '-', s_char(2)
139			WRITE(, ) s_pos(1), '-', f_pos(1), '\|', f_pos(2), '-', s_pos(2)
140			ENDIF
141			!--
142			IF (nbops > nbops_max-1) THEN
143			CALL histerr(3, 'decoop', 'Expression too complex', str, ' ')
144			ENDIF
145			!--
146			IF (check) WRITE(, ) 'decoop : --', nbops, ' ', str(1:leng)
147			!---
148			!-- Start the analysis of the syntax. 3 types of constructs
149			!-- are recognized. They are scanned sequentialy
150			!---
151			IF (nbsep == 1) THEN
152			IF (check) WRITE(, ) 'decoop : Only one operation'
153			IF (INDEX(ops, f_char(1)) > 0) THEN
154			!------ Type : scal+X
155			IF (f_char(1) == '-' .OR. f_char(1) == '/') THEN
156			opp_str = f_char(1)//'I'
157			ELSE
158			opp_str = f_char(1)
159			ENDIF
160			scal_str = str(s_pos(1)+1:f_pos(1)-1)
161			str = 'X'
162			ELSE IF (INDEX(ops, f_char(2)) > 0) THEN
163			!------ Type : X+scal
164			opp_str = f_char(2)
165			scal_str = str(f_pos(2)+1:s_pos(2)-1)
166			str = 'X'
167			ELSE
168			CALL histerr(3, 'decoop', &
169			& 'Unknown operations of type X+scal', f_char(1), pstr)
170			ENDIF
171			ELSE
172			IF (check) WRITE(, ) 'decoop : More complex operation'
173			IF ( f_char(1) == '(' .AND. f_char(2) == ')' ) THEN
174			!------ Type : sin(X)
175			opp_str = str(s_pos(1)+1:f_pos(1)-1)
176			scal_str = '?'
177			str = str(1:s_pos(1))//'X'//str(f_pos(2)+1:leng)
178			ELSE IF ( (f_char(1) == '(' .AND. f_char(2) == ', ')&
179			& .OR.(f_char(1) == ', ' .AND. f_char(2) == ')')) THEN
180			!------ Type : max(X, scal) or max(scal, X)
181			IF (f_char(1) == '(' .AND. s_char(2) == ')') THEN
182			!-------- Type : max(X, scal)
183			opp_str = str(f_pos(1)-3:f_pos(1)-1)
184			scal_str = str(f_pos(2)+1:s_pos(2)-1)
185			str = str(1:f_pos(1)-4)//'X'//str(s_pos(2)+1:leng)
186			ELSE IF (f_char(1) == ', ' .AND. s_char(1) == '(') THEN
187			!-------- Type : max(scal, X)
188			opp_str = str(s_pos(1)-3:s_pos(1)-1)
189			scal_str = str(s_pos(1)+1:f_pos(1)-1)
190			str = str(1:s_pos(1)-4)//'X'//str(f_pos(2)+1:leng)
191			ELSE
192			CALL histerr(3, 'decoop', 'Syntax error 1', str, ' ')
193			ENDIF
194			ELSE
195			prio = (f_char(2) == '*').OR.(f_char(2) == '^')
196			IF ( (INDEX(ops, f_char(1)) > 0) &
197			& .AND.(xpos-f_pos(1) == 1).AND.(.NOT.prio) ) THEN
198			!-------- Type : ... scal+X ...
199			IF (f_char(1) == '-' .OR. f_char(1) == '/') THEN
200			opp_str = f_char(1)//'I'
201			ELSE
202			opp_str = f_char(1)
203			ENDIF
204			scal_str = str(s_pos(1)+1:f_pos(1)-1)
205			str = str(1:s_pos(1))//'X'//str(f_pos(1)+2:leng)
206			ELSE IF ( (INDEX(ops, f_char(2)) > 0) &
207			& .AND.(f_pos(2)-xpos == 1) ) THEN
208			!-------- Type : ... X+scal ...
209			opp_str = f_char(2)
210			scal_str = str(f_pos(2)+1:s_pos(2)-1)
211			str = str(1:f_pos(2)-2)//'X'//str(s_pos(2):leng)
212			ELSE
213			CALL histerr(3, 'decoop', 'Syntax error 2', str, ' ')
214			ENDIF
215			ENDIF
216			ENDIF
217			!---
218			IF (check) WRITE(, ) 'decoop : Finished syntax, str = ', TRIM(str)
219			!---
220			!-- Now that the different components of the operation are identified
221			!-- we transform them into what is going to be used in the program
222			!---
223			IF (INDEX(scal_str, '?') > 0) THEN
224			IF (INDEX(funcs, opp_str(1:LEN_TRIM(opp_str))) > 0) THEN
225			opps(nbops) = opp_str(1:LEN_TRIM(opp_str))
226			scal(nbops) = missing_val
227			ELSE
228			CALL histerr(3, 'decoop', &
229			& 'Unknown function', opp_str(1:LEN_TRIM(opp_str)), ' ')
230			ENDIF
231			ELSE
232			leng = LEN_TRIM(opp_str)
233			IF (INDEX(mima, opp_str(1:leng)) > 0) THEN
234			opps(nbops) = 'fu'//opp_str(1:leng)
235			ELSE
236			IF (INDEX(opp_str(1:leng), '+') > 0) THEN
237			opps(nbops) = 'add'
238			ELSE IF (INDEX(opp_str(1:leng), '-I') > 0) THEN
239			opps(nbops) = 'subi'
240			ELSE IF (INDEX(opp_str(1:leng), '-') > 0) THEN
241			opps(nbops) = 'sub'
242			ELSE IF (INDEX(opp_str(1:leng), '*') > 0) THEN
243			opps(nbops) = 'mult'
244			ELSE IF (INDEX(opp_str(1:leng), '/') > 0) THEN
245			opps(nbops) = 'div'
246			ELSE IF (INDEX(opp_str(1:leng), '/I') > 0) THEN
247			opps(nbops) = 'divi'
248			ELSE IF (INDEX(opp_str(1:leng), '^') > 0) THEN
249			opps(nbops) = 'power'
250			ELSE
251			CALL histerr(3, 'decoop', &
252			& 'Unknown operation', opp_str(1:leng), ' ')
253			ENDIF
254			ENDIF
255			!-----
256			leng = LEN_TRIM(scal_str)
257			ppos = INDEX(scal_str, '.')
258			epos = INDEX(scal_str, 'e')
259			IF (epos == 0) epos = INDEX(scal_str, 'E')
260			!-----
261			!---- Try to catch a few errors
262			!-----
263			IF (INDEX(ops, scal_str) > 0) THEN
264			CALL histerr(3, 'decoop', &
265			& 'Strange scalar you have here ', scal_str, pstr)
266			ENDIF
267			IF (epos > 0) THEN
268			WRITE(tl, '(I3.3)') leng
269			WRITE(dl, '(I3.3)') epos-ppos-1
270			fmt='(e'//tl//'.'//dl//')'
271			READ(scal_str, fmt) scal(nbops)
272			ELSE IF (ppos > 0) THEN
273			WRITE(tl, '(I3.3)') leng
274			WRITE(dl, '(I3.3)') leng-ppos
275			fmt='(f'//tl//'.'//dl//')'
276			READ(scal_str, fmt) scal(nbops)
277			ELSE
278			WRITE(tl, '(I3.3)') leng
279			fmt = '(I'//tl//')'
280			READ(scal_str, fmt) int_tmp
281			scal(nbops) = REAL(int_tmp)
282			ENDIF
283			ENDIF
284			IF (check) WRITE(, ) 'decoop : Finished interpretation'
285			CALL findsep(str, nbsep, f_char, f_pos, s_char, s_pos)
286			ENDDO
287			!--------------------
288			END SUBROUTINE decoop
289
290			!************************************************
291
292			SUBROUTINE findsep (str, nbsep, f_char, f_pos, s_char, s_pos)
293			!- Subroutine finds all separators in a given string
294			!- It returns the following information about str :
295			!- f_char : The first separation character
296			!- (1 for before and 2 for after)
297			!- f_pos : The position of the first separator
298			!- s_char : The second separation character
299			!- (1 for before and 2 for after)
300			!- s_pos : The position of the second separator
301	guez	32	USE errioipsl, ONLY : histerr
302
303	guez	30	CHARACTER(LEN=80) :: str
304			INTEGER :: nbsep
305			CHARACTER(LEN=1), DIMENSION(2) :: f_char, s_char
306			INTEGER, DIMENSION(2) :: f_pos, s_pos
307
308			CHARACTER(LEN=70) :: str_tmp
309			LOGICAL :: f_found, s_found
310			INTEGER :: ind, xpos, leng, i
311
312			LOGICAL :: check = .FALSE.
313			!---------------------------------------------------------------------
314			IF (check) WRITE(, ) 'findsep : call cleanstr: ', TRIM(str)
315
316			CALL cleanstr(str)
317
318			IF (check) WRITE(, ) 'findsep : out of cleanstr: ', TRIM(str)
319
320			xpos = INDEX(str, 'X')
321			leng = LEN_TRIM(str)
322
323			f_pos(1:2) = (/ 0, leng+1 /)
324			f_char(1:2) = (/ '?', '?' /)
325			s_pos(1:2) = (/ 0, leng+1 /)
326			s_char(1:2) = (/ '?', '?' /)
327
328			nbsep = 0
329
330			f_found = .FALSE.
331			s_found = .FALSE.
332			IF (xpos > 1) THEN
333			DO i=xpos-1, 1, -1
334			ind = INDEX(seps, str(i:i))
335			IF (ind > 0) THEN
336			IF (.NOT.f_found) THEN
337			f_char(1) = str(i:i)
338			f_pos(1) = i
339			nbsep = nbsep+1
340			f_found = .TRUE.
341			ELSE IF (.NOT.s_found) THEN
342			s_char(1) = str(i:i)
343			s_pos(1) = i
344			nbsep = nbsep+1
345			s_found = .TRUE.
346			ENDIF
347			ENDIF
348			ENDDO
349			ENDIF
350
351			f_found = .FALSE.
352			s_found = .FALSE.
353			IF (xpos < leng) THEN
354			DO i=xpos+1, leng
355			ind = INDEX(seps, str(i:i))
356			IF (ind > 0) THEN
357			IF (.NOT.f_found) THEN
358			f_char(2) = str(i:i)
359			f_pos(2) = i
360			nbsep = nbsep+1
361			f_found = .TRUE.
362			ELSE IF (.NOT.s_found) THEN
363			s_char(2) = str(i:i)
364			s_pos(2) = i
365			nbsep = nbsep+1
366			s_found = .TRUE.
367			ENDIF
368			ENDIF
369			ENDDO
370			ENDIF
371
372			IF (nbsep > 4) THEN
373			WRITE(str_tmp, '("number :", I3)') nbsep
374			CALL histerr(3, 'findsep', &
375			& 'How can I find that many separators', str_tmp, str)
376			ENDIF
377
378			IF (check) WRITE(, ) 'Finished findsep : ', nbsep, leng
379			!---------------------
380			END SUBROUTINE findsep
381
382			!************************************************
383
384			SUBROUTINE cleanstr(str)
385			!- We clean up the string by taking out the extra () and puting
386			!- everything in lower case except for the X describing the variable
387	guez	32	use strlowercase_m, only: strlowercase
388
389	guez	30	CHARACTER(LEN=80) :: str
390
391			INTEGER :: ind, leng, ic, it
392			LOGICAL :: check = .FALSE.
393			!---------------------------------------------------------------------
394			leng = LEN_TRIM(str)
395			CALL strlowercase(str)
396
397			ind = INDEX(str, 'x')
398			IF (check) THEN
399			WRITE (, ) 'cleanstr 1.0 : ind = ', ind, &
400			& ' str = ', str(1:leng), '---'
401			ENDIF
402
403			! If the character before the x is not a letter then we can assume
404			! that it is the variable and promote it to a capital letter
405
406			DO WHILE (ind > 0)
407			ic = 0
408			IF (ind > 1) ic = IACHAR(str(ind-1:ind-1))
409			IF (ic < 97 .OR. ic > 122) THEN
410			str(ind:ind) = 'X'
411			ENDIF
412			it = INDEX(str(ind+1:leng), 'x')
413			IF (it > 0) THEN
414			ind = ind+it
415			ELSE
416			ind = it
417			ENDIF
418			ENDDO
419
420			IF (check) WRITE (, ) 'cleanstr 2.0 : str = ', str(1:leng), '---'
421
422			IF ( str(1:1) == '(' .AND. str(leng:leng) == ')' ) THEN
423			str = str(2:leng-1)
424			ENDIF
425
426			IF (check) WRITE (, ) 'cleanstr 3.0 : str = ', str(1:leng), '---'
427
428			leng = LEN_TRIM(str)
429			ind = INDEX(str, '((X))')
430			IF (ind > 0) THEN
431			str=str(1:ind-1)//'(X)'//str(ind+5:leng)//' '
432			ENDIF
433
434			IF (check) WRITE (, ) 'cleanstr 4.0 : str = ', str(1:leng), '---'
435
436			leng = LEN_TRIM(str)
437			ind = INDEX(str, '(X)')
438			IF (ind > 0 .AND. ind+3 < leng) THEN
439			IF ( (INDEX(seps, str(ind-1:ind-1)) > 0) &
440			& .AND. (INDEX(seps, str(ind+3:ind+3)) > 0) ) THEN
441			str=str(1:ind-1)//'X'//str(ind+3:leng)//' '
442			ENDIF
443			ENDIF
444
445			IF (check) WRITE (, ) 'cleanstr 5.0 : str = ', str(1:leng), '---'
446
447			leng = LEN_TRIM(str)
448			ind = INDEX(str(1:leng), ' ')
449			DO WHILE (ind > 0)
450			str=str(1:ind-1)//str(ind+1:leng)//' '
451			leng = LEN_TRIM(str)
452			ind = INDEX(str(1:leng), ' ')
453			ENDDO
454
455			IF (check) WRITE (, ) 'cleanstr 6.0 : str = ', str(1:leng), '---'
456			!----------------------
457			END SUBROUTINE cleanstr
458
459			!************************************************
460
461			SUBROUTINE moycum (opp, np, px, py, pwx)
462			!- Does time operations
463	guez	32	USE errioipsl, ONLY : histerr
464
465	guez	30	CHARACTER(LEN=7) :: opp
466			INTEGER :: np
467			REAL, DIMENSION(:) :: px, py
468			INTEGER :: pwx
469			!---------------------------------------------------------------------
470			IF (pwx /= 0) THEN
471			IF (opp == 'ave') THEN
472			px(1:np)=(px(1:np)*pwx+py(1:np))/REAL(pwx+1)
473			ELSE IF (opp == 't_sum') THEN
474			px(1:np)=px(1:np)+py(1:np)
475			ELSE IF ( (opp == 'l_min').OR.(opp == 't_min') ) THEN
476			px(1:np)=MIN(px(1:np), py(1:np))
477			ELSE IF ( (opp == 'l_max').OR.(opp == 't_max') ) THEN
478			px(1:np)=MAX(px(1:np), py(1:np))
479			ELSE
480			CALL histerr(3, "moycum", 'Unknown time operation', opp, ' ')
481			ENDIF
482			ELSE
483			IF (opp == 'l_min') THEN
484			px(1:np)=MIN(px(1:np), py(1:np))
485			ELSE IF (opp == 'l_max') THEN
486			px(1:np)=MAX(px(1:np), py(1:np))
487			ELSE
488			px(1:np)=py(1:np)
489			ENDIF
490			ENDIF
491			!--------------------
492			END SUBROUTINE moycum
493
494			!************************************************
495
496			SUBROUTINE trans_buff (ox, sx, oy, sy, oz, sz, xsz, ysz, zsz, v3d, sl, v1d)
497			!- This subroutine extracts from the full 3D variable the slab of
498			!- data that will be used later. Perhaps there are hardware routines
499			!- for this task on some computers. This routine will be obsolete in
500			!- a F90 environnement
501
502			!- INPUT
503			!- ox : Origin of slab of data in X
504			!- sx : Size of slab in X
505			!- oy : Origin of slab of data in Y
506			!- sy : Size of slab in Y
507			!- oz : Origin of slab of data in Z
508			!- sz : Size of slab in Z
509			!- xsz, ysz, zsz : 3 sizes of full variable v3d
510			!- v3d : The full 3D variable
511			!- sl : size of variable v1d
512			!- v1d : The 1D variable containing the slab
513
514			INTEGER :: ox, sx, oy, sy, oz, sz
515			INTEGER :: xsz, ysz, zsz
516			INTEGER :: sl
517			REAL :: v3d(xsz, ysz, zsz)
518			REAL :: v1d(sl)
519
520			!---------------------------------------------------------------------
521
522			V1d(:sxsysz) = reshape(V3d(ox:ox-1+sx, oy:oy-1+sy, oz:oz-1+sz), &
523			SHAPE=(/sxsysz/))
524
525			END SUBROUTINE trans_buff
526
527			END MODULE mathelp