MODULE mathelp
!-
!$Id: mathelp.f90 1927 2012-11-22 13:54:21Z dsolyga $
!-
! This software is governed by the CeCILL license
! See IOIPSL/IOIPSL_License_CeCILL.txt
!---------------------------------------------------------------------
  USE errioipsl,ONLY : ipslerr
  USE stringop
!-
  PRIVATE
  PUBLIC :: mathop,moycum,moycum_index,buildop
!-
  INTERFACE mathop
    MODULE PROCEDURE mathop_r11,mathop_r21,mathop_r31
  END INTERFACE
!-
!- Variables used to detect and identify the operations
!-
  CHARACTER(LEN=80),SAVE :: &
 &  seps='( ) , + - / * ^', ops = '+ - * / ^', 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'
  CHARACTER(LEN=120),SAVE :: &
 &  indexfu = 'gather, scatter, fill, coll, undef, only'
!---------------------------------------------------------------------
CONTAINS
!===
SUBROUTINE buildop (c_str,ex_topps,topp,fill_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
!-
!- c_str    : String containing the operations
!- ex_toops : Time operations that can be expected within the string
!- fill_val :
!-
!- OUTPUT
!-
!- topp   : Time operation
!- opps   :
!- scal   :
!- nbops  :
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  CHARACTER(LEN=*),INTENT(IN) :: c_str,ex_topps
  CHARACTER(LEN=*),INTENT(OUT) :: topp
  CHARACTER(LEN=*),DIMENSION(:),INTENT(OUT) :: opps
  REAL,INTENT(IN) :: fill_val
  REAL,DIMENSION(:),INTENT(OUT) :: scal
  INTEGER,INTENT(OUT) :: nbops
!-
  CHARACTER(LEN=LEN(c_str)) :: str,new_str
  INTEGER :: leng,ind_opb,ind_clb
!-
  LOGICAL :: check = .FALSE.
!---------------------------------------------------------------------
  IF (check) WRITE(*,*) 'buildop : Some preliminary cleaning'
!-
  str = c_str
  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,fill_val,opps,scal,nbops)
      ELSE
        CALL ipslerr(3,'buildop', &
 &        'time operation does not exist',str(1:ind_opb-1),' ')
      ENDIF
    ELSE
      CALL ipslerr(3,'buildop', &
 &      'some long operation 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,fill_val,opps,scal,nbops)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  CHARACTER(LEN=*),INTENT(IN) :: pstr
  REAL,INTENT(IN) :: fill_val
  CHARACTER(LEN=*),DIMENSION(:),INTENT(OUT) :: opps
  REAL,DIMENSION(:),INTENT(OUT) :: scal
  INTEGER,INTENT(OUT) :: nbops
!-
  CHARACTER(LEN=1),DIMENSION(2) :: f_char,s_char
  INTEGER,DIMENSION(2) :: f_pos,s_pos
  CHARACTER(LEN=20) :: opp_str,scal_str
  CHARACTER(LEN=LEN(pstr)) :: str
  INTEGER :: nbsep,nbops_max,xpos,leng,ppos,epos,int_tmp
  CHARACTER(LEN=3) :: tl,dl
  CHARACTER(LEN=10) :: fmt
!-
  LOGICAL :: check = .FALSE.,prio
!---------------------------------------------------------------------
  IF (check) WRITE(*,'(2A)') ' decoop : Incoming string : ',pstr
!-
  str = pstr; nbops = 0;
!-
  CALL findsep (str,nbsep,f_char,f_pos,s_char,s_pos)
  IF (check) WRITE(*,*) 'decoop : Out of findsep',nbsep
!-
  nbops_max = min(SIZE(opps),SIZE(scal))
!-
  DO WHILE (nbsep > 0)
    IF (nbops >= nbops_max) THEN
      CALL ipslerr(3,'decoop','Expression too complex',TRIM(str),' ')
    ENDIF
!--
    xpos = INDEX(str,'X')
    leng = LEN_TRIM(str)
    nbops = nbops+1
!--
    IF (check) THEN
      WRITE(*,*) 'decoop : str   -> ',TRIM(str)
      WRITE(*,*) 'decoop : nbops -> ',nbops
      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
!---
!-- 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 ipslerr(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 ipslerr(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 ipslerr(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) = fill_val
      ELSE
        CALL ipslerr(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 ipslerr(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 ipslerr(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
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  CHARACTER(LEN=*),INTENT(INOUT) :: str
  INTEGER :: nbsep
  CHARACTER(LEN=1),DIMENSION(2) :: f_char,s_char
  INTEGER,DIMENSION(2) :: f_pos,s_pos
!-
  CHARACTER(LEN=10) :: 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 ipslerr(3,'findsep', &
 &    'How can I find that many separators',str_tmp,TRIM(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
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  CHARACTER(LEN=*),INTENT(INOUT) :: 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 mathop_r11 &
 &  (fun,nb,work_in,miss_val,nb_index,nindex,scal,nb_max,work_out)
!---------------------------------------------------------------------
!- This subroutines gives an interface to the various operation
!- which are allowed. The interface is general enough to allow its use
!- for other cases.
!-
!- INPUT
!-
!- fun      : function to be applied to the vector of data
!- nb       : Length of input vector
!- work_in  : Input vector of data (REAL)
!- miss_val : The value of the missing data flag (it has to be a
!-            maximum value, in f90 : huge( a real ))
!- nb_index : Length of index vector
!- nindex   : Vector of indices
!- scal     : A scalar value for vector/scalar operations
!- nb_max   : maximum length of output vector
!-
!- OUTPUT
!-
!- nb_max   : Actual length of output variable
!- work_out : Output vector after the operation was applied
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  CHARACTER(LEN=7) :: fun
  INTEGER :: nb,nb_max,nb_index
  INTEGER :: nindex(nb_index)
  REAL :: work_in(nb),scal,miss_val
  REAL :: work_out(nb_max)
!-
  INTEGER :: ierr
!-
  INTRINSIC SIN,COS,TAN,ASIN,ACOS,ATAN,EXP,LOG,SQRT,ABS
!---------------------------------------------------------------------
  ierr = 0
!-
  IF (scal >= miss_val-1.) THEN
    IF (INDEX(indexfu,fun(1:LEN_TRIM(fun))) == 0) THEN
      SELECT CASE (fun)
      CASE('sin')
        ierr = ma_sin_r11(nb,work_in,nb_max,work_out)
      CASE('cos')
        ierr = ma_cos_r11(nb,work_in,nb_max,work_out)
      CASE('tan')
        ierr = ma_tan_r11(nb,work_in,nb_max,work_out)
      CASE('asin')
        ierr = ma_asin_r11(nb,work_in,nb_max,work_out)
      CASE('acos')
        ierr = ma_acos_r11(nb,work_in,nb_max,work_out)
      CASE('atan')
        ierr = ma_atan_r11(nb,work_in,nb_max,work_out)
      CASE('exp')
        ierr = ma_exp_r11(nb,work_in,nb_max,work_out)
      CASE('log')
        ierr = ma_log_r11(nb,work_in,nb_max,work_out)
      CASE('sqrt')
        ierr = ma_sqrt_r11(nb,work_in,nb_max,work_out)
      CASE('chs')
        ierr = ma_chs_r11(nb,work_in,nb_max,work_out)
      CASE('abs')
        ierr = ma_abs_r11(nb,work_in,nb_max,work_out)
      CASE('cels')
        ierr = ma_cels_r11(nb,work_in,nb_max,work_out)
      CASE('kelv')
        ierr = ma_kelv_r11(nb,work_in,nb_max,work_out)
      CASE('deg')
        ierr = ma_deg_r11(nb,work_in,nb_max,work_out)
      CASE('rad')
        ierr = ma_rad_r11(nb,work_in,nb_max,work_out)
      CASE('ident')
        ierr = ma_ident_r11(nb,work_in,nb_max,work_out)
      CASE DEFAULT
        CALL ipslerr(3,"mathop", &
 &        'scalar variable undefined and no indexing', &
 &        'but still unknown function',fun)
      END SELECT
      IF (ierr > 0) THEN
        CALL ipslerr(3,"mathop", &
 &        'Error while executing a simple function',fun,' ')
      ENDIF
    ELSE
      SELECT CASE (fun)
      CASE('gather')
        ierr = ma_fugath_r11(nb,work_in,nb_index,nindex, &
&                            miss_val,nb_max,work_out)
      CASE('scatter')
        IF (nb_index > nb) THEN
          work_out(1:nb_max) = miss_val
          ierr=1
        ELSE
          ierr = ma_fuscat_r11(nb,work_in,nb_index,nindex, &
&                              miss_val,nb_max,work_out)
        ENDIF
      CASE('coll')
        ierr = ma_fucoll_r11(nb,work_in,nb_index,nindex, &
&                            miss_val,nb_max,work_out)
      CASE('fill')
        ierr = ma_fufill_r11(nb,work_in,nb_index,nindex, &
&                            miss_val,nb_max,work_out)
      CASE('undef')
        ierr = ma_fuundef_r11(nb,work_in,nb_index,nindex, &
&                            miss_val,nb_max,work_out)
      CASE('only')
        ierr = ma_fuonly_r11(nb,work_in,nb_index,nindex, &
&                            miss_val,nb_max,work_out)
      CASE DEFAULT
        CALL ipslerr(3,"mathop", &
 &        'scalar variable undefined and indexing',&
 &        'was requested but with unknown function',fun)
      END SELECT
      IF (ierr > 0) THEN
        CALL ipslerr(3,"mathop_r11", &
 &        'Error while executing an indexing function',fun,' ')
      ENDIF
    ENDIF
  ELSE
    SELECT CASE (fun)
    CASE('fumin')
      ierr = ma_fumin_r11(nb,work_in,scal,nb_max,work_out)
    CASE('fumax')
      ierr = ma_fumax_r11(nb,work_in,scal,nb_max,work_out)
    CASE('add')
      ierr = ma_add_r11(nb,work_in,scal,nb_max,work_out)
    CASE('subi')
      ierr = ma_subi_r11(nb,work_in,scal,nb_max,work_out)
    CASE('sub')
      ierr = ma_sub_r11(nb,work_in,scal,nb_max,work_out)
    CASE('mult')
      ierr = ma_mult_r11(nb,work_in,scal,nb_max,work_out)
    CASE('div')
      ierr = ma_div_r11(nb,work_in,scal,nb_max,work_out)
    CASE('divi')
      ierr = ma_divi_r11(nb,work_in,scal,nb_max,work_out)
    CASE('power')
      ierr = ma_power_r11(nb,work_in,scal,nb_max,work_out)
    CASE DEFAULT
      CALL ipslerr(3,"mathop", &
 &      'Unknown operation with a scalar',fun,' ')
    END SELECT
    IF (ierr > 0) THEN
      CALL ipslerr(3,"mathop", &
 &      'Error while executing a scalar function',fun,' ')
    ENDIF
  ENDIF
!------------------------
END SUBROUTINE mathop_r11
!-
!=== FUNCTIONS (only one argument)
!-
INTEGER FUNCTION ma_sin_r11(nb,x,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb,nbo,i
  REAL :: x(nb),y(nbo)
!---------------------------------------------------------------------
  DO i=1,nb
    y(i) = SIN(x(i))
  ENDDO
!-
  nbo = nb
  ma_sin_r11 = 0
!----------------------
END FUNCTION ma_sin_r11
!===
INTEGER FUNCTION ma_cos_r11(nb,x,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb,nbo,i
  REAL :: x(nb),y(nbo)
!---------------------------------------------------------------------
  DO i=1,nb
    y(i) = COS(x(i))
  ENDDO
!-
  nbo = nb
  ma_cos_r11 = 0
!----------------------
END FUNCTION ma_cos_r11
!===
INTEGER FUNCTION ma_tan_r11(nb,x,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb,nbo,i
  REAL :: x(nb),y(nbo)
!---------------------------------------------------------------------
  DO i=1,nb
    y(i) = TAN(x(i))
  ENDDO
!-
  nbo = nb
  ma_tan_r11 = 0
!----------------------
END FUNCTION ma_tan_r11
!===
INTEGER FUNCTION ma_asin_r11(nb,x,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb,nbo,i
  REAL :: x(nb),y(nbo)
!---------------------------------------------------------------------
  DO i=1,nb
    y(i) = ASIN(x(i))
  ENDDO
!-
  nbo = nb
  ma_asin_r11 = 0
!-----------------------
END FUNCTION ma_asin_r11
!===
INTEGER FUNCTION ma_acos_r11(nb,x,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb,nbo,i
  REAL :: x(nb),y(nbo)
!---------------------------------------------------------------------
  DO i=1,nb
    y(i) = ACOS(x(i))
  ENDDO
!-
  nbo = nb
  ma_acos_r11 = 0
!-----------------------
END FUNCTION ma_acos_r11
!===
INTEGER FUNCTION ma_atan_r11(nb,x,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb,nbo,i
  REAL :: x(nb),y(nbo)
!---------------------------------------------------------------------
  DO i=1,nb
    y(i) = ATAN(x(i))
  ENDDO
!-
  nbo = nb
  ma_atan_r11 = 0
!-----------------------
END FUNCTION ma_atan_r11
!===
INTEGER FUNCTION ma_exp_r11(nb,x,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb,nbo,i
  REAL :: x(nb),y(nbo)
!---------------------------------------------------------------------
  DO i=1,nb
    y(i) = EXP(x(i))
  ENDDO
!-
  nbo = nb
  ma_exp_r11 = 0
!----------------------
END FUNCTION ma_exp_r11
!===
INTEGER FUNCTION ma_log_r11(nb,x,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb,nbo,i
  REAL :: x(nb),y(nbo)
!---------------------------------------------------------------------
  DO i=1,nb
    y(i) = log(x(i))
  ENDDO
!-
  nbo = nb
  ma_log_r11 = 0
!----------------------
END FUNCTION ma_log_r11
!===
INTEGER FUNCTION ma_sqrt_r11(nb,x,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb,nbo,i
  REAL :: x(nb),y(nbo)
!---------------------------------------------------------------------
  DO i=1,nb
    y(i) = SQRT(x(i))
  ENDDO
!-
  nbo = nb
  ma_sqrt_r11 = 0
!-----------------------
END FUNCTION ma_sqrt_r11
!===
INTEGER FUNCTION ma_abs_r11(nb,x,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb,nbo,i
  REAL :: x(nb),y(nbo)
!---------------------------------------------------------------------
  DO i=1,nb
    y(i) = ABS(x(i))
  ENDDO
!-
  nbo = nb
  ma_abs_r11 = 0
!----------------------
END FUNCTION ma_abs_r11
!===
INTEGER FUNCTION ma_chs_r11(nb,x,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb,nbo,i
  REAL :: x(nb),y(nbo)
!---------------------------------------------------------------------
  DO i=1,nb
    y(i) = x(i)*(-1.)
  ENDDO
!-
  nbo = nb
  ma_chs_r11 = 0
!----------------------
END FUNCTION ma_chs_r11
!===
INTEGER FUNCTION ma_cels_r11(nb,x,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb,nbo,i
  REAL :: x(nb),y(nbo)
!---------------------------------------------------------------------
  DO i=1,nb
    y(i) = x(i)-273.15
  ENDDO
!-
  nbo = nb
  ma_cels_r11 = 0
!-----------------------
END FUNCTION ma_cels_r11
!===
INTEGER FUNCTION ma_kelv_r11(nb,x,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb,nbo,i
  REAL :: x(nb),y(nbo)
!---------------------------------------------------------------------
  DO i=1,nb
    y(i) = x(i)+273.15
  ENDDO
!-
  nbo = nb
  ma_kelv_r11 = 0
!-----------------------
END FUNCTION ma_kelv_r11
!===
INTEGER FUNCTION ma_deg_r11(nb,x,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb,nbo,i
  REAL :: x(nb),y(nbo)
!---------------------------------------------------------------------
  DO i=1,nb
    y(i) = x(i)*57.29577951
  ENDDO
!-
  nbo = nb
  ma_deg_r11 = 0
!-----------------------
END FUNCTION ma_deg_r11
!===
INTEGER FUNCTION ma_rad_r11(nb,x,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb,nbo,i
  REAL :: x(nb),y(nbo)
!---------------------------------------------------------------------
  DO i=1,nb
    y(i) = x(i)*0.01745329252
  ENDDO
!-
  nbo = nb
  ma_rad_r11 = 0
!----------------------
END FUNCTION ma_rad_r11
!===
INTEGER FUNCTION ma_ident_r11(nb,x,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb,nbo,i
  REAL :: x(nb),y(nbo)
!---------------------------------------------------------------------
  DO i=1,nb
    y(i) = x(i)
  ENDDO
!-
  nbo = nb
  ma_ident_r11 = 0
!------------------------
END FUNCTION ma_ident_r11
!-
!=== OPERATIONS (two argument)
!-
INTEGER FUNCTION ma_add_r11(nb,x,s,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb,nbo
  REAL :: x(nb),s,y(nbo)
!-
  INTEGER :: i
!---------------------------------------------------------------------
  DO i=1,nb
    y(i) = x(i)+s
  ENDDO
!-
  nbo = nb
  ma_add_r11 = 0
!-----------------------
  END FUNCTION ma_add_r11
!===
INTEGER FUNCTION ma_sub_r11(nb,x,s,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb,nbo
  REAL :: x(nb),s,y(nbo)
!-
  INTEGER :: i
!---------------------------------------------------------------------
  DO i=1,nb
    y(i) = x(i)-s
  ENDDO
!-
  nbo = nb
  ma_sub_r11 = 0
!----------------------
END FUNCTION ma_sub_r11
!===
INTEGER FUNCTION ma_subi_r11(nb,x,s,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb,nbo
  REAL :: x(nb),s,y(nbo)
!-
  INTEGER :: i
!---------------------------------------------------------------------
  DO i=1,nb
    y(i) =  s-x(i)
  ENDDO
!-
  nbo = nb
  ma_subi_r11 = 0
!-----------------------
END FUNCTION ma_subi_r11
!===
INTEGER FUNCTION ma_mult_r11(nb,x,s,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb,nbo
  REAL :: x(nb),s,y(nbo)
!-
  INTEGER :: i
!---------------------------------------------------------------------
  DO i=1,nb
    y(i) = x(i)*s
  ENDDO
!-
  nbo = nb
  ma_mult_r11 = 0
!-----------------------
END FUNCTION ma_mult_r11
!===
INTEGER FUNCTION ma_div_r11(nb,x,s,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb,nbo
  REAL :: x(nb),s,y(nbo)
!-
  INTEGER :: i
!---------------------------------------------------------------------
  DO i=1,nb
    y(i) = x(i)/s
  ENDDO
!-
  nbo = nb
  ma_div_r11 = 0
!-----------------------
  END FUNCTION ma_div_r11
!===
INTEGER FUNCTION ma_divi_r11(nb,x,s,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb,nbo
  REAL :: x(nb),s,y(nbo)
!-
  INTEGER :: i
!---------------------------------------------------------------------
  DO i=1,nb
    y(i) = s/x(i)
  ENDDO
!-
  nbo = nb
  ma_divi_r11 = 0
!-----------------------
END FUNCTION ma_divi_r11
!===
INTEGER FUNCTION ma_power_r11(nb,x,s,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb,nbo
  REAL :: x(nb),s,y(nbo)
!-
  INTEGER :: i
!---------------------------------------------------------------------
  DO i=1,nb
    y(i) = x(i)**s
  ENDDO
!-
  nbo = nb
  ma_power_r11 = 0
!-----------------------
END FUNCTION ma_power_r11
!===
INTEGER FUNCTION ma_fumin_r11(nb,x,s,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb,nbo
  REAL :: x(nb),s,y(nbo)
!-
  INTEGER :: i
!---------------------------------------------------------------------
  DO i=1,nb
    y(i) = MIN(x(i),s)
  ENDDO
!-
  nbo = nb
  ma_fumin_r11 = 0
!------------------------
END FUNCTION ma_fumin_r11
!===
INTEGER FUNCTION ma_fumax_r11(nb,x,s,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb,nbo
  REAL :: x(nb),s,y(nbo)
!-
  INTEGER :: i
!---------------------------------------------------------------------
  DO i=1,nb
    y(i) = MAX(x(i),s)
  ENDDO
!-
  nbo = nb
  ma_fumax_r11 = 0
!------------------------
END FUNCTION ma_fumax_r11
!===
INTEGER FUNCTION ma_fuscat_r11(nb,x,nbi,ind,miss_val,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb,nbo,nbi
  INTEGER :: ind(nbi)
  REAL :: x(nb),miss_val,y(nbo)
!-
  INTEGER :: i,ii,ipos
!---------------------------------------------------------------------
  ma_fuscat_r11 = 0
!-
  y(1:nbo) = miss_val
!-
  IF (nbi <= nb) THEN
    ipos = 0
    DO i=1,nbi
      IF (ind(i) <= nbo .AND. ind(i) > 0) THEN
        ipos = ipos+1
        y(ind(i)) = x(ipos)
      ELSE
        IF (ind(i) > nbo) ma_fuscat_r11  = ma_fuscat_r11+1
      ENDIF
    ENDDO
!-- Repeat the data if needed
    IF (MINVAL(ind) < 0) THEN
      DO i=1,nbi
        IF (ind(i) <= 0) THEN
          DO ii=1,ABS(ind(i))-1
            IF (ind(i+1)+ii <= nbo) THEN
              y(ind(i+1)+ii) = y(ind(i+1))
            ELSE
              ma_fuscat_r11  = ma_fuscat_r11+1
            ENDIF
          ENDDO
        ENDIF
      ENDDO
    ENDIF
  ELSE
    ma_fuscat_r11  = 1
  ENDIF
!-------------------------
END FUNCTION ma_fuscat_r11
!===
INTEGER FUNCTION ma_fugath_r11(nb,x,nbi,ind,miss_val,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb,nbo,nbi
  INTEGER :: ind(nbi)
  REAL :: x(nb),miss_val,y(nbo)
!-
  INTEGER :: i,ipos
!---------------------------------------------------------------------
  IF (nbi <= nbo) THEN
    ma_fugath_r11 = 0
    y(1:nbo) = miss_val
    ipos = 0
    DO i=1,nbi
      IF (ipos+1 <= nbo) THEN
        IF (ind(i) > 0) THEN
          ipos = ipos+1
          y(ipos) = x(ind(i))
        ENDIF
      ELSE
        IF (ipos+1 > nbo) ma_fugath_r11  = ma_fugath_r11+1
      ENDIF
    ENDDO
  ELSE
    ma_fugath_r11 = 1
  ENDIF
!-
  nbo = ipos
!-------------------------
END FUNCTION ma_fugath_r11
!===
INTEGER FUNCTION ma_fufill_r11(nb,x,nbi,ind,miss_val,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb,nbo,nbi
  INTEGER :: ind(nbi)
  REAL :: x(nb),miss_val,y(nbo)
!-
  INTEGER :: i,ii,ipos
!---------------------------------------------------------------------
  ma_fufill_r11 = 0
!-
  IF (nbi <= nb) THEN
    ipos = 0
    DO i=1,nbi
      IF (ind(i) <= nbo .AND. ind(i) > 0) THEN
        ipos = ipos+1
        y(ind(i)) = x(ipos)
      ELSE
        IF (ind(i) > nbo) ma_fufill_r11  = ma_fufill_r11+1
      ENDIF
    ENDDO
!-- Repeat the data if needed
    IF (MINVAL(ind) < 0) THEN
      DO i=1,nbi
        IF (ind(i) <= 0) THEN
          DO ii=1,ABS(ind(i))-1
            IF (ind(i+1)+ii <= nbo) THEN
              y(ind(i+1)+ii) = y(ind(i+1))
            ELSE
              ma_fufill_r11  = ma_fufill_r11+1
            ENDIF
          ENDDO
        ENDIF
      ENDDO
    ENDIF
  ELSE
    ma_fufill_r11  = 1
  ENDIF
!-------------------------
END FUNCTION ma_fufill_r11
!===
INTEGER FUNCTION ma_fucoll_r11(nb,x,nbi,ind,miss_val,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb,nbo,nbi
  INTEGER :: ind(nbi)
  REAL :: x(nb),miss_val,y(nbo)
!-
  INTEGER :: i,ipos
!---------------------------------------------------------------------
  IF (nbi <= nbo) THEN
    ma_fucoll_r11 = 0
    ipos = 0
    DO i=1,nbi
      IF (ipos+1 <= nbo) THEN
        IF (ind(i) > 0) THEN
          ipos = ipos+1
          y(ipos) = x(ind(i))
        ENDIF
      ELSE
        IF (ipos+1 > nbo) ma_fucoll_r11  = ma_fucoll_r11+1
      ENDIF
    ENDDO
  ELSE
    ma_fucoll_r11 = 1
  ENDIF
!-
  nbo = ipos
!-------------------------
END FUNCTION ma_fucoll_r11
!===
INTEGER FUNCTION ma_fuundef_r11(nb,x,nbi,ind,miss_val,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb,nbo,nbi
  INTEGER :: ind(nbi)
  REAL :: x(nb),miss_val,y(nbo)
!-
  INTEGER :: i
!---------------------------------------------------------------------
  IF (nbi <= nbo .AND. nbo == nb) THEN
    ma_fuundef_r11 = 0
    DO i=1,nbo
      y(i) = x(i)
    ENDDO
    DO i=1,nbi
      IF (ind(i) <= nbo .AND. ind(i) > 0) THEN
        y(ind(i)) =  miss_val
      ELSE
        IF (ind(i) > nbo) ma_fuundef_r11  = ma_fuundef_r11+1
      ENDIF
    ENDDO
  ELSE
    ma_fuundef_r11 = 1
  ENDIF
!--------------------------
END FUNCTION ma_fuundef_r11
!===
INTEGER FUNCTION ma_fuonly_r11(nb,x,nbi,ind,miss_val,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb,nbo,nbi
  INTEGER :: ind(nbi)
  REAL :: x(nb),miss_val,y(nbo)
!-
  INTEGER :: i
!---------------------------------------------------------------------
  IF (     (nbi <= nbo).AND.(nbo == nb) &
 &    .AND.ALL(ind(1:nbi) <= nbo) ) THEN
    ma_fuonly_r11 = 0
    y(1:nbo) = miss_val
    DO i=1,nbi
      IF (ind(i) > 0) THEN
        y(ind(i)) =  x(ind(i))
      ENDIF
    ENDDO
  ELSE
    ma_fuonly_r11 = 1
  ENDIF
!-------------------------
END FUNCTION ma_fuonly_r11
!===
!===
SUBROUTINE mathop_r21 &
 &  (fun,nb,work_in,miss_val,nb_index,nindex,scal,nb_max,work_out)
!---------------------------------------------------------------------
!- This subroutines gives an interface to the various operations
!- which are allowed. The interface is general enough to allow its use
!- for other cases.
!-
!- INPUT
!-
!- fun      : function to be applied to the vector of data
!- nb       : Length of input vector
!- work_in  : Input vector of data (REAL)
!- miss_val : The value of the missing data flag (it has to be a
!-            maximum value, in f90 : huge( a real ))
!- nb_index : Length of index vector
!- nindex   : Vector of indices
!- scal     : A scalar value for vector/scalar operations
!- nb_max   : maximum length of output vector
!-
!- OUTPUT
!-
!- nb_max   : Actual length of output variable
!- work_out : Output vector after the operation was applied
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  CHARACTER(LEN=7) :: fun
  INTEGER :: nb(2),nb_max,nb_index
  INTEGER :: nindex(nb_index)
  REAL :: work_in(nb(1),nb(2)),scal,miss_val
  REAL :: work_out(nb_max)
!-
  INTEGER :: ierr
!-
  INTRINSIC SIN,COS,TAN,ASIN,ACOS,ATAN,EXP,LOG,SQRT,ABS
!---------------------------------------------------------------------
  ierr = 0
!-
  IF (scal >= miss_val-1.) THEN
    IF (INDEX(indexfu,fun(1:LEN_TRIM(fun))) == 0) THEN
      SELECT CASE (fun)
      CASE('sin')
        ierr = ma_sin_r21(nb,work_in,nb_max,work_out)
      CASE('cos')
        ierr = ma_cos_r21(nb,work_in,nb_max,work_out)
      CASE('tan')
        ierr = ma_tan_r21(nb,work_in,nb_max,work_out)
      CASE('asin')
        ierr = ma_asin_r21(nb,work_in,nb_max,work_out)
      CASE('acos')
        ierr = ma_acos_r21(nb,work_in,nb_max,work_out)
      CASE('atan')
        ierr = ma_atan_r21(nb,work_in,nb_max,work_out)
      CASE('exp')
        ierr = ma_exp_r21(nb,work_in,nb_max,work_out)
      CASE('log')
        ierr = ma_log_r21(nb,work_in,nb_max,work_out)
      CASE('sqrt')
        ierr = ma_sqrt_r21(nb,work_in,nb_max,work_out)
      CASE('chs')
        ierr = ma_chs_r21(nb,work_in,nb_max,work_out)
      CASE('abs')
        ierr = ma_abs_r21(nb,work_in,nb_max,work_out)
      CASE('cels')
        ierr = ma_cels_r21(nb,work_in,nb_max,work_out)
      CASE('kelv')
        ierr = ma_kelv_r21(nb,work_in,nb_max,work_out)
      CASE('deg')
        ierr = ma_deg_r21(nb,work_in,nb_max,work_out)
      CASE('rad')
        ierr = ma_rad_r21(nb,work_in,nb_max,work_out)
      CASE('ident')
        ierr = ma_ident_r21(nb,work_in,nb_max,work_out)
      CASE DEFAULT
        CALL ipslerr(3,"mathop", &
 &        'scalar variable undefined and no indexing', &
 &        'but still unknown function',fun)
      END SELECT
      IF (ierr > 0) THEN
        CALL ipslerr(3,"mathop", &
 &        'Error while executing a simple function',fun,' ')
      ENDIF
    ELSE
      SELECT CASE (fun)
      CASE('gather')
        ierr = ma_fugath_r21(nb,work_in,nb_index,nindex, &
 &                            miss_val,nb_max,work_out)
      CASE('scatter')
        IF (nb_index > (nb(1)*nb(2)) ) THEN
          work_out(1:nb_max) = miss_val
          ierr=1
        ELSE
          ierr = ma_fuscat_r21(nb,work_in,nb_index,nindex, &
 &                             miss_val,nb_max,work_out)
        ENDIF
      CASE('coll')
        ierr = ma_fucoll_r21(nb,work_in,nb_index,nindex, &
 &                           miss_val,nb_max,work_out)
      CASE('fill')
        ierr = ma_fufill_r21(nb,work_in,nb_index,nindex, &
 &                           miss_val,nb_max,work_out)
      CASE('undef')
        ierr = ma_fuundef_r21(nb,work_in,nb_index,nindex, &
 &                           miss_val,nb_max,work_out)
      CASE('only')
        ierr = ma_fuonly_r21(nb,work_in,nb_index,nindex, &
 &                           miss_val,nb_max,work_out)
      CASE DEFAULT
        CALL ipslerr(3,"mathop", &
 &        'scalar variable undefined and indexing', &
 &        'was requested but with unknown function',fun)
      END SELECT
      IF (ierr > 0) THEN
        CALL ipslerr(3,"mathop_r21", &
 &        'Error while executing an indexing function',fun,' ')
      ENDIF
    ENDIF
  ELSE
    SELECT CASE (fun)
    CASE('fumin')
      ierr = ma_fumin_r21(nb,work_in,scal,nb_max,work_out)
    CASE('fumax')
      ierr = ma_fumax_r21(nb,work_in,scal,nb_max,work_out)
    CASE('add')
      ierr = ma_add_r21(nb,work_in,scal,nb_max,work_out)
    CASE('subi')
      ierr = ma_subi_r21(nb,work_in,scal,nb_max,work_out)
    CASE('sub')
      ierr = ma_sub_r21(nb,work_in,scal,nb_max,work_out)
    CASE('mult')
      ierr = ma_mult_r21(nb,work_in,scal,nb_max,work_out)
    CASE('div')
      ierr = ma_div_r21(nb,work_in,scal,nb_max,work_out)
    CASE('divi')
      ierr = ma_divi_r21(nb,work_in,scal,nb_max,work_out)
    CASE('power')
      ierr = ma_power_r21(nb,work_in,scal,nb_max,work_out)
    CASE DEFAULT
      CALL ipslerr(3,"mathop", &
 &      'Unknown operation with a scalar',fun,' ')
    END SELECT
    IF (ierr > 0) THEN
      CALL ipslerr(3,"mathop", &
 &      'Error while executing a scalar function',fun,' ')
    ENDIF
  ENDIF
!------------------------
END SUBROUTINE mathop_r21
!-
!=== FUNCTIONS (only one argument)
!-
INTEGER FUNCTION ma_sin_r21(nb,x,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(2),nbo,i,j,ij
  REAL :: x(nb(1),nb(2)),y(nbo)
!---------------------------------------------------------------------
  ij = 0
  DO j=1,nb(2)
    DO i=1,nb(1)
      ij = ij+1
      y(ij) = SIN(x(i,j))
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)
  ma_sin_r21 = 0
!----------------------
END FUNCTION ma_sin_r21
!===
INTEGER FUNCTION ma_cos_r21(nb,x,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(2),nbo,i,j,ij
  REAL :: x(nb(1),nb(2)),y(nbo)
!---------------------------------------------------------------------
  ij = 0
  DO j=1,nb(2)
    DO i=1,nb(1)
      ij = ij+1
      y(ij) = COS(x(i,j))
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)
  ma_cos_r21 = 0
!----------------------
END FUNCTION ma_cos_r21
!===
INTEGER FUNCTION ma_tan_r21(nb,x,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(2),nbo,i,j,ij
  REAL :: x(nb(1),nb(2)),y(nbo)
!---------------------------------------------------------------------
  ij = 0
  DO j=1,nb(2)
    DO i=1,nb(1)
      ij = ij+1
      y(ij) = TAN(x(i,j))
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)
  ma_tan_r21 = 0
!----------------------
END FUNCTION ma_tan_r21
!===
  INTEGER FUNCTION ma_asin_r21(nb,x,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(2),nbo,i,j,ij
  REAL :: x(nb(1),nb(2)),y(nbo)
!---------------------------------------------------------------------
  ij = 0
  DO j=1,nb(2)
    DO i=1,nb(1)
      ij = ij+1
      y(ij) = ASIN(x(i,j))
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)
  ma_asin_r21 = 0
!-----------------------
END FUNCTION ma_asin_r21
!===
INTEGER FUNCTION ma_acos_r21(nb,x,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(2),nbo,i,j,ij
  REAL :: x(nb(1),nb(2)),y(nbo)
!---------------------------------------------------------------------
  ij = 0
  DO j=1,nb(2)
    DO i=1,nb(1)
      ij = ij+1
      y(ij) = ACOS(x(i,j))
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)
  ma_acos_r21 = 0
!-----------------------
END FUNCTION ma_acos_r21
!===
INTEGER FUNCTION ma_atan_r21(nb,x,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(2),nbo,i,j,ij
  REAL :: x(nb(1),nb(2)),y(nbo)
!---------------------------------------------------------------------
  ij = 0
  DO j=1,nb(2)
    DO i=1,nb(1)
      ij = ij+1
      y(ij) = ATAN(x(i,j))
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)
  ma_atan_r21 = 0
!-----------------------
END FUNCTION ma_atan_r21
!===
INTEGER FUNCTION ma_exp_r21(nb,x,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(2),nbo,i,j,ij
  REAL :: x(nb(1),nb(2)),y(nbo)
!---------------------------------------------------------------------
  ij = 0
  DO j=1,nb(2)
    DO i=1,nb(1)
      ij = ij+1
      y(ij) = EXP(x(i,j))
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)
  ma_exp_r21 = 0
!----------------------
END FUNCTION ma_exp_r21
!===
INTEGER FUNCTION ma_log_r21(nb,x,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(2),nbo,i,j,ij
  REAL :: x(nb(1),nb(2)),y(nbo)
!---------------------------------------------------------------------
  ij = 0
  DO j=1,nb(2)
    DO i=1,nb(1)
      ij = ij+1
      y(ij) = LOG(x(i,j))
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)
  ma_log_r21 = 0
!----------------------
END FUNCTION ma_log_r21
!===
INTEGER FUNCTION ma_sqrt_r21(nb,x,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(2),nbo,i,j,ij
  REAL :: x(nb(1),nb(2)),y(nbo)
!---------------------------------------------------------------------
  ij = 0
  DO j=1,nb(2)
    DO i=1,nb(1)
      ij = ij+1
      y(ij) = SQRT(x(i,j))
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)
  ma_sqrt_r21 = 0
!-----------------------
END FUNCTION ma_sqrt_r21
!===
INTEGER FUNCTION ma_abs_r21(nb,x,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(2),nbo,i,j,ij
  REAL :: x(nb(1),nb(2)),y(nbo)
!---------------------------------------------------------------------
  ij = 0
  DO j=1,nb(2)
    DO i=1,nb(1)
      ij = ij+1
      y(ij) = ABS(x(i,j))
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)
  ma_abs_r21 = 0
!----------------------
END FUNCTION ma_abs_r21
!===
INTEGER FUNCTION ma_chs_r21(nb,x,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(2),nbo,i,j,ij
  REAL :: x(nb(1),nb(2)),y(nbo)
!---------------------------------------------------------------------
  ij = 0
  DO j=1,nb(2)
    DO i=1,nb(1)
      ij = ij+1
      y(ij) = x(i,j)*(-1.)
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)
  ma_chs_r21 = 0
!----------------------
END FUNCTION ma_chs_r21
!===
INTEGER FUNCTION ma_cels_r21(nb,x,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(2),nbo,i,j,ij
  REAL :: x(nb(1),nb(2)),y(nbo)
!---------------------------------------------------------------------
  ij = 0
  DO j=1,nb(2)
    DO i=1,nb(1)
      ij = ij+1
      y(ij) = x(i,j)-273.15
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)
  ma_cels_r21 = 0
!-----------------------
END FUNCTION ma_cels_r21
!===
INTEGER FUNCTION ma_kelv_r21(nb,x,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(2),nbo,i,j,ij
  REAL :: x(nb(1),nb(2)),y(nbo)
!---------------------------------------------------------------------
  ij = 0
  DO j=1,nb(2)
    DO i=1,nb(1)
      ij = ij+1
      y(ij) = x(i,j)+273.15
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)
  ma_kelv_r21 = 0
!-----------------------
END FUNCTION ma_kelv_r21
!===
INTEGER FUNCTION ma_deg_r21(nb,x,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(2),nbo,i,j,ij
  REAL :: x(nb(1),nb(2)),y(nbo)
!---------------------------------------------------------------------
  ij = 0
  DO j=1,nb(2)
    DO i=1,nb(1)
      ij = ij+1
      y(ij) = x(i,j)*57.29577951
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)
  ma_deg_r21 = 0
!----------------------
END FUNCTION ma_deg_r21
!===
INTEGER FUNCTION ma_rad_r21(nb,x,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(2),nbo,i,j,ij
  REAL :: x(nb(1),nb(2)),y(nbo)
!---------------------------------------------------------------------
  ij = 0
  DO j=1,nb(2)
    DO i=1,nb(1)
      ij = ij+1
      y(ij) = x(i,j)*0.01745329252
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)
  ma_rad_r21 = 0
!----------------------
END FUNCTION ma_rad_r21
!===
INTEGER FUNCTION ma_ident_r21(nb,x,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(2),nbo,i,j,ij
  REAL :: x(nb(1),nb(2)),y(nbo)
!---------------------------------------------------------------------
  ij = 0
  DO j=1,nb(2)
    DO i=1,nb(1)
      ij = ij+1
      y(ij) = x(i,j)
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)
  ma_ident_r21 = 0
!------------------------
END FUNCTION ma_ident_r21
!-
!=== OPERATIONS (two argument)
!-
INTEGER FUNCTION ma_add_r21(nb,x,s,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(2),nbo
  REAL :: x(nb(1),nb(2)),s,y(nbo)
!-
  INTEGER :: i,j,ij
!---------------------------------------------------------------------
  ij = 0
  DO j=1,nb(2)
    DO i=1,nb(1)
      ij = ij+1
      y(ij) = x(i,j)+s
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)
  ma_add_r21 = 0
!----------------------
END FUNCTION ma_add_r21
!===
INTEGER FUNCTION ma_sub_r21(nb,x,s,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(2),nbo
  REAL :: x(nb(1),nb(2)),s,y(nbo)
!-
  INTEGER :: i,j,ij
!---------------------------------------------------------------------
  ij = 0
  DO j=1,nb(2)
    DO i=1,nb(1)
      ij = ij+1
      y(ij) = x(i,j)-s
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)
  ma_sub_r21 = 0
!----------------------
END FUNCTION ma_sub_r21
!===
INTEGER FUNCTION ma_subi_r21(nb,x,s,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(2),nbo
  REAL :: x(nb(1),nb(2)),s,y(nbo)
!-
  INTEGER :: i,j,ij
!---------------------------------------------------------------------
  ij = 0
  DO j=1,nb(2)
    DO i=1,nb(1)
      ij = ij+1
      y(ij) =  s-x(i,j)
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)
  ma_subi_r21 = 0
!-----------------------
END FUNCTION ma_subi_r21
!===
INTEGER FUNCTION ma_mult_r21(nb,x,s,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(2),nbo
  REAL :: x(nb(1),nb(2)),s,y(nbo)
!-
  INTEGER :: i,j,ij
!---------------------------------------------------------------------
  ij = 0
  DO j=1,nb(2)
    DO i=1,nb(1)
      ij = ij+1
      y(ij) = x(i,j)*s
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)
  ma_mult_r21 = 0
!-----------------------
END FUNCTION ma_mult_r21
!===
INTEGER FUNCTION ma_div_r21(nb,x,s,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(2),nbo
  REAL :: x(nb(1),nb(2)),s,y(nbo)
!-
  INTEGER :: i,j,ij
!---------------------------------------------------------------------
  ij = 0
  DO j=1,nb(2)
    DO i=1,nb(1)
      ij = ij+1
      y(ij) = x(i,j)/s
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)
  ma_div_r21 = 0
!----------------------
END FUNCTION ma_div_r21
!===
INTEGER FUNCTION ma_divi_r21(nb,x,s,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(2),nbo
  REAL :: x(nb(1),nb(2)),s,y(nbo)
!-
  INTEGER :: i,j,ij
!---------------------------------------------------------------------
  ij = 0
  DO j=1,nb(2)
    DO i=1,nb(1)
      ij = ij+1
      y(ij) = s/x(i,j)
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)
  ma_divi_r21 = 0
!-----------------------
END FUNCTION ma_divi_r21
!===
INTEGER FUNCTION ma_power_r21(nb,x,s,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(2),nbo
  REAL :: x(nb(1),nb(2)),s,y(nbo)
!-
  INTEGER :: i,j,ij
!---------------------------------------------------------------------
  ij = 0
  DO j=1,nb(2)
    DO i=1,nb(1)
      ij = ij+1
      y(ij) = x(i,j) ** s
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)
  ma_power_r21 = 0
!------------------------
END FUNCTION ma_power_r21
!===
INTEGER FUNCTION ma_fumin_r21(nb,x,s,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(2),nbo
  REAL :: x(nb(1),nb(2)),s,y(nbo)
!-
  INTEGER :: i,j,ij
!---------------------------------------------------------------------
  ij = 0
  DO j=1,nb(2)
    DO i=1,nb(1)
      ij = ij+1
      y(ij) = MIN(x(i,j),s)
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)
  ma_fumin_r21 = 0
!------------------------
END FUNCTION ma_fumin_r21
!===
INTEGER FUNCTION ma_fumax_r21(nb,x,s,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(2),nbo
  REAL :: x(nb(1),nb(2)),s,y(nbo)
!-
  INTEGER :: i,j,ij
!---------------------------------------------------------------------
  ij = 0
  DO j=1,nb(2)
    DO i=1,nb(1)
      ij = ij+1
      y(ij) = MAX(x(i,j),s)
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)
  ma_fumax_r21 = 0
!------------------------
END FUNCTION ma_fumax_r21
!===
INTEGER FUNCTION ma_fuscat_r21(nb,x,nbi,ind,miss_val,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(2),nbo,nbi
  INTEGER :: ind(nbi)
  REAL :: x(nb(1),nb(2)),miss_val,y(nbo)
!-
  INTEGER :: i,j,ij,ii,ipos
!---------------------------------------------------------------------
  ma_fuscat_r21 = 0
!-
  y(1:nbo) = miss_val
!-
  IF (nbi <= nb(1)*nb(2)) THEN
    ipos = 0
    DO ij=1,nbi
       IF (ind(ij) <= nbo .AND. ind(ij) > 0) THEN
         ipos = ipos+1
         j = ((ipos-1)/nb(1))+1
         i = (ipos-(j-1)*nb(1))
         y(ind(ij)) = x(i,j)
       ELSE
         IF (ind(ij) > nbo) ma_fuscat_r21  = ma_fuscat_r21+1
       ENDIF
    ENDDO
!-- Repeat the data if needed
    IF (MINVAL(ind) < 0) THEN
      DO i=1,nbi
        IF (ind(i) <= 0) THEN
          DO ii=1,ABS(ind(i))-1
            IF (ind(i+1)+ii <= nbo) THEN
              y(ind(i+1)+ii) = y(ind(i+1))
            ELSE
              ma_fuscat_r21  = ma_fuscat_r21+1
            ENDIF
          ENDDO
        ENDIF
      ENDDO
    ENDIF
  ELSE
    ma_fuscat_r21  = 1
  ENDIF
!-------------------------
END FUNCTION ma_fuscat_r21
!===
INTEGER FUNCTION ma_fugath_r21(nb,x,nbi,ind,miss_val,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(2),nbo,nbi
  INTEGER :: ind(nbi)
  REAL :: x(nb(1),nb(2)),miss_val,y(nbo)
!-
  INTEGER :: i,j,ij,ipos
!---------------------------------------------------------------------
  IF (nbi <= nbo) THEN
    ma_fugath_r21 = 0
    y(1:nbo) = miss_val
    ipos = 0
    DO ij=1,nbi
      IF (ipos+1 <= nbo) THEN
        IF (ind(ij) > 0) THEN
          j = ((ind(ij)-1)/nb(1))+1
          i = (ind(ij)-(j-1)*nb(1))
          ipos = ipos+1
          y(ipos) = x(i,j)
        ENDIF
      ELSE
        IF (ipos+1 > nbo) ma_fugath_r21  = ma_fugath_r21+1
      ENDIF
    ENDDO
  ELSE
    ma_fugath_r21 = 1
  ENDIF
  nbo = ipos
!-------------------------
END FUNCTION ma_fugath_r21
!===
INTEGER FUNCTION ma_fufill_r21(nb,x,nbi,ind,miss_val,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(2),nbo,nbi
  INTEGER :: ind(nbi)
  REAL :: x(nb(1),nb(2)),miss_val,y(nbo)
!-
  INTEGER :: i,j,ij,ii,ipos
!---------------------------------------------------------------------
  ma_fufill_r21 = 0
!-
  IF (nbi <= nb(1)*nb(2)) THEN
    ipos = 0
    DO ij=1,nbi
      IF (ind(ij) <= nbo .AND. ind(ij) > 0) THEN
        ipos = ipos+1
        j = ((ipos-1)/nb(1))+1
        i = (ipos-(j-1)*nb(1))
        y(ind(ij)) = x(i,j)
      ELSE
        IF (ind(ij) > nbo) ma_fufill_r21  = ma_fufill_r21+1
      ENDIF
    ENDDO
!-- Repeat the data if needed
    IF (MINVAL(ind) < 0) THEN
      DO i=1,nbi
        IF (ind(i) <= 0) THEN
          DO ii=1,ABS(ind(i))-1
            IF (ind(i+1)+ii <= nbo) THEN
              y(ind(i+1)+ii) = y(ind(i+1))
            ELSE
              ma_fufill_r21  = ma_fufill_r21+1
            ENDIF
          ENDDO
        ENDIF
      ENDDO
    ENDIF
  ELSE
    ma_fufill_r21  = 1
  ENDIF
!-------------------------
END FUNCTION ma_fufill_r21
!===
INTEGER FUNCTION ma_fucoll_r21(nb,x,nbi,ind,miss_val,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(2),nbo,nbi
  INTEGER :: ind(nbi)
  REAL :: x(nb(1),nb(2)),miss_val,y(nbo)
!-
  INTEGER :: i,j,ij,ipos
!---------------------------------------------------------------------
  IF (nbi <= nbo) THEN
    ma_fucoll_r21 = 0
    ipos = 0
    DO ij=1,nbi
      IF (ipos+1 <= nbo) THEN
        IF (ind(ij) > 0) THEN
          j = ((ind(ij)-1)/nb(1))+1
          i = (ind(ij)-(j-1)*nb(1))
          ipos = ipos+1
          y(ipos) = x(i,j)
        ENDIF
      ELSE
        IF (ipos+1 > nbo) ma_fucoll_r21  = ma_fucoll_r21+1
      ENDIF
    ENDDO
  ELSE
    ma_fucoll_r21 = 1
  ENDIF
  nbo = ipos
!-------------------------
END FUNCTION ma_fucoll_r21
!===
INTEGER FUNCTION ma_fuundef_r21(nb,x,nbi,ind,miss_val,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(2),nbo,nbi
  INTEGER :: ind(nbi)
  REAL :: x(nb(1),nb(2)),miss_val,y(nbo)
!-
  INTEGER :: i,j,ij
!---------------------------------------------------------------------
  IF (nbi <= nbo .AND. nbo == nb(1)*nb(2)) THEN
    ma_fuundef_r21 = 0
    DO ij=1,nbo
      j = ((ij-1)/nb(1))+1
      i = (ij-(j-1)*nb(1))
      y(ij) = x(i,j)
    ENDDO
    DO i=1,nbi
      IF (ind(i) <= nbo .AND. ind(i) > 0) THEN
        y(ind(i)) =  miss_val
      ELSE
        IF (ind(i) > nbo) ma_fuundef_r21 = ma_fuundef_r21+1
      ENDIF
    ENDDO
  ELSE
    ma_fuundef_r21 = 1
  ENDIF
!--------------------------
END FUNCTION ma_fuundef_r21
!===
INTEGER FUNCTION ma_fuonly_r21(nb,x,nbi,ind,miss_val,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(2),nbo,nbi
  INTEGER :: ind(nbi)
  REAL :: x(nb(1),nb(2)),miss_val,y(nbo)
!-
  INTEGER :: i,j,ij
!---------------------------------------------------------------------
  IF (     (nbi <= nbo).AND.(nbo == nb(1)*nb(2)) &
 &    .AND.ALL(ind(1:nbi) <= nbo) ) THEN
    ma_fuonly_r21 = 0
    y(1:nbo) = miss_val
    DO ij=1,nbi
      IF (ind(ij) > 0) THEN
        j = ((ind(ij)-1)/nb(1))+1
        i = (ind(ij)-(j-1)*nb(1))
        y(ind(ij)) =  x(i,j)
      ENDIF
    ENDDO
  ELSE
    ma_fuonly_r21 = 1
  ENDIF
!-------------------------
END FUNCTION ma_fuonly_r21
!===
!===
SUBROUTINE mathop_r31 &
 &  (fun,nb,work_in,miss_val,nb_index,nindex,scal,nb_max,work_out)
!---------------------------------------------------------------------
!- This subroutines gives an interface to the various operations
!- which are allowed. The interface is general enough to allow its use
!- for other cases.
!-
!- INPUT
!-
!- fun      : function to be applied to the vector of data
!- nb       : Length of input vector
!- work_in  : Input vector of data (REAL)
!- miss_val : The value of the missing data flag (it has to be a
!-            maximum value, in f90 : huge( a real ))
!- nb_index : Length of index vector
!- nindex   : Vector of indices
!- scal     : A scalar value for vector/scalar operations
!- nb_max   : maximum length of output vector
!-
!- OUTPUT
!-
!- nb_max   : Actual length of output variable
!- work_out : Output vector after the operation was applied
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  CHARACTER(LEN=7) :: fun
  INTEGER :: nb(3),nb_max,nb_index
  INTEGER :: nindex(nb_index)
  REAL :: work_in(nb(1),nb(2),nb(3)),scal,miss_val
  REAL :: work_out(nb_max)
!-
  INTEGER :: ierr
!-
  INTRINSIC SIN,COS,TAN,ASIN,ACOS,ATAN,EXP,LOG,SQRT,ABS
!---------------------------------------------------------------------
  ierr = 0
!-
  IF (scal >= miss_val-1.) THEN
    IF (INDEX(indexfu,fun(1:LEN_TRIM(fun))) == 0) THEN
      SELECT CASE (fun)
      CASE('sin')
        ierr = ma_sin_r31(nb,work_in,nb_max,work_out)
      CASE('cos')
        ierr = ma_cos_r31(nb,work_in,nb_max,work_out)
      CASE('tan')
        ierr = ma_tan_r31(nb,work_in,nb_max,work_out)
      CASE('asin')
        ierr = ma_asin_r31(nb,work_in,nb_max,work_out)
      CASE('acos')
        ierr = ma_acos_r31(nb,work_in,nb_max,work_out)
      CASE('atan')
        ierr = ma_atan_r31(nb,work_in,nb_max,work_out)
      CASE('exp')
        ierr = ma_exp_r31(nb,work_in,nb_max,work_out)
      CASE('log')
        ierr = ma_log_r31(nb,work_in,nb_max,work_out)
      CASE('sqrt')
        ierr = ma_sqrt_r31(nb,work_in,nb_max,work_out)
      CASE('chs')
        ierr = ma_chs_r31(nb,work_in,nb_max,work_out)
      CASE('abs')
        ierr = ma_abs_r31(nb,work_in,nb_max,work_out)
      CASE('cels')
        ierr = ma_cels_r31(nb,work_in,nb_max,work_out)
      CASE('kelv')
        ierr = ma_kelv_r31(nb,work_in,nb_max,work_out)
      CASE('deg')
        ierr = ma_deg_r31(nb,work_in,nb_max,work_out)
      CASE('rad')
        ierr = ma_rad_r31(nb,work_in,nb_max,work_out)
      CASE('ident')
        ierr = ma_ident_r31(nb,work_in,nb_max,work_out)
      CASE DEFAULT
        CALL ipslerr(3,"mathop", &
 &        'scalar variable undefined and no indexing', &
 &        'but still unknown function',fun)
      END SELECT
      IF (ierr > 0) THEN
         CALL ipslerr(3,"mathop", &
 &        'Error while executing a simple function',fun,' ')
      ENDIF
    ELSE
      SELECT CASE (fun)
      CASE('gather')
        ierr = ma_fugath_r31(nb,work_in,nb_index,nindex, &
 &                           miss_val,nb_max,work_out)
      CASE('scatter')
        IF (nb_index > (nb(1)*nb(2)*nb(3))) THEN
          work_out(1:nb_max) = miss_val
          ierr=1
        ELSE
          ierr = ma_fuscat_r31(nb,work_in,nb_index,nindex, &
 &                             miss_val,nb_max,work_out)
        ENDIF
      CASE('coll')
        ierr = ma_fucoll_r31(nb,work_in,nb_index,nindex, &
 &                           miss_val,nb_max,work_out)
      CASE('fill')
        ierr = ma_fufill_r31(nb,work_in,nb_index,nindex, &
 &                           miss_val,nb_max,work_out)
      CASE('undef')
        ierr = ma_fuundef_r31(nb,work_in,nb_index,nindex, &
 &                           miss_val,nb_max,work_out)
      CASE('only')
        ierr = ma_fuonly_r31(nb,work_in,nb_index,nindex, &
 &                           miss_val,nb_max,work_out)
      CASE DEFAULT
        CALL ipslerr(3,"mathop", &
 &        'scalar variable undefined and indexing', &
 &        'was requested but with unknown function',fun)
      END SELECT
      IF (ierr > 0) THEN
        CALL ipslerr(3,"mathop_r31", &
 &        'Error while executing an indexing function',fun,' ')
      ENDIF
    ENDIF
  ELSE
    SELECT CASE (fun)
    CASE('fumin')
      ierr = ma_fumin_r31(nb,work_in,scal,nb_max,work_out)
    CASE('fumax')
      ierr = ma_fumax_r31(nb,work_in,scal,nb_max,work_out)
    CASE('add')
      ierr = ma_add_r31(nb,work_in,scal,nb_max,work_out)
    CASE('subi')
      ierr = ma_subi_r31(nb,work_in,scal,nb_max,work_out)
    CASE('sub')
      ierr = ma_sub_r31(nb,work_in,scal,nb_max,work_out)
    CASE('mult')
      ierr = ma_mult_r31(nb,work_in,scal,nb_max,work_out)
    CASE('div')
      ierr = ma_div_r31(nb,work_in,scal,nb_max,work_out)
    CASE('divi')
      ierr = ma_divi_r31(nb,work_in,scal,nb_max,work_out)
    CASE('power')
      ierr = ma_power_r31(nb,work_in,scal,nb_max,work_out)
    CASE DEFAULT
      CALL ipslerr(3,"mathop", &
 &      'Unknown operation with a scalar',fun,' ')
    END SELECT
    IF (ierr > 0) THEN
      CALL ipslerr(3,"mathop", &
 &      'Error while executing a scalar function',fun,' ')
    ENDIF
  ENDIF
!------------------------
END SUBROUTINE mathop_r31
!-
!=== FUNCTIONS (only one argument)
!-
INTEGER FUNCTION ma_sin_r31(nb,x,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(3),nbo,i,j,k,ij
  REAL :: x(nb(1),nb(2),nb(3)),y(nbo)
!---------------------------------------------------------------------
  ij = 0
  DO k=1,nb(3)
    DO j=1,nb(2)
      DO i=1,nb(1)
        ij = ij+1
        y(ij) = SIN(x(i,j,k))
      ENDDO
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)*nb(3)
  ma_sin_r31 = 0
!----------------------
END FUNCTION ma_sin_r31
!===
INTEGER FUNCTION ma_cos_r31(nb,x,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(3),nbo,i,j,k,ij
  REAL :: x(nb(1),nb(2),nb(3)),y(nbo)
!---------------------------------------------------------------------
  ij = 0
  DO k=1,nb(3)
    DO j=1,nb(2)
      DO i=1,nb(1)
        ij = ij+1
        y(ij) = COS(x(i,j,k))
      ENDDO
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)*nb(3)
  ma_cos_r31 = 0
!----------------------
END FUNCTION ma_cos_r31
!===
INTEGER FUNCTION ma_tan_r31(nb,x,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(3),nbo,i,j,k,ij
  REAL :: x(nb(1),nb(2),nb(3)),y(nbo)
!---------------------------------------------------------------------
  ij = 0
  DO k=1,nb(3)
    DO j=1,nb(2)
      DO i=1,nb(1)
        ij = ij+1
        y(ij) = TAN(x(i,j,k))
      ENDDO
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)*nb(3)
  ma_tan_r31 = 0
!----------------------
END FUNCTION ma_tan_r31
!===
INTEGER FUNCTION ma_asin_r31(nb,x,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(3),nbo,i,j,k,ij
  REAL :: x(nb(1),nb(2),nb(3)),y(nbo)
!---------------------------------------------------------------------
  ij = 0
  DO k=1,nb(3)
    DO j=1,nb(2)
      DO i=1,nb(1)
        ij = ij+1
        y(ij) = ASIN(x(i,j,k))
      ENDDO
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)*nb(3)
  ma_asin_r31 = 0
!-----------------------
END FUNCTION ma_asin_r31
!===
INTEGER FUNCTION ma_acos_r31(nb,x,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(3),nbo,i,j,k,ij
  REAL :: x(nb(1),nb(2),nb(3)),y(nbo)
!---------------------------------------------------------------------
  ij = 0
  DO k=1,nb(3)
    DO j=1,nb(2)
      DO i=1,nb(1)
        ij = ij+1
        y(ij) = ACOS(x(i,j,k))
      ENDDO
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)*nb(3)
  ma_acos_r31 = 0
!-----------------------
END FUNCTION ma_acos_r31
!===
INTEGER FUNCTION ma_atan_r31(nb,x,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(3),nbo,i,j,k,ij
  REAL :: x(nb(1),nb(2),nb(3)),y(nbo)
!---------------------------------------------------------------------
  ij = 0
  DO k=1,nb(3)
    DO j=1,nb(2)
      DO i=1,nb(1)
        ij = ij+1
        y(ij) = ATAN(x(i,j,k))
      ENDDO
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)*nb(3)
  ma_atan_r31 = 0
!-----------------------
  END FUNCTION ma_atan_r31
!===
INTEGER FUNCTION ma_exp_r31(nb,x,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(3),nbo,i,j,k,ij
  REAL :: x(nb(1),nb(2),nb(3)),y(nbo)
!---------------------------------------------------------------------
  ij = 0
  DO k=1,nb(3)
    DO j=1,nb(2)
      DO i=1,nb(1)
        ij = ij+1
        y(ij) = EXP(x(i,j,k))
      ENDDO
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)*nb(3)
  ma_exp_r31 = 0
!----------------------
END FUNCTION ma_exp_r31
!===
INTEGER FUNCTION ma_log_r31(nb,x,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(3),nbo,i,j,k,ij
  REAL :: x(nb(1),nb(2),nb(3)),y(nbo)
!---------------------------------------------------------------------
  ij = 0
  DO k=1,nb(3)
    DO j=1,nb(2)
      DO i=1,nb(1)
        ij = ij+1
        y(ij) = LOG(x(i,j,k))
      ENDDO
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)*nb(3)
  ma_log_r31 = 0
!----------------------
END FUNCTION ma_log_r31
!===
INTEGER FUNCTION ma_sqrt_r31(nb,x,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(3),nbo,i,j,k,ij
  REAL :: x(nb(1),nb(2),nb(3)),y(nbo)
!---------------------------------------------------------------------
  ij = 0
  DO k=1,nb(3)
    DO j=1,nb(2)
      DO i=1,nb(1)
        ij = ij+1
        y(ij) = SQRT(x(i,j,k))
      ENDDO
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)*nb(3)
  ma_sqrt_r31 = 0
!-----------------------
END FUNCTION ma_sqrt_r31
!===
INTEGER FUNCTION ma_abs_r31(nb,x,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(3),nbo,i,j,k,ij
  REAL :: x(nb(1),nb(2),nb(3)),y(nbo)
!---------------------------------------------------------------------
  ij = 0
  DO k=1,nb(3)
    DO j=1,nb(2)
      DO i=1,nb(1)
        ij = ij+1
        y(ij) = ABS(x(i,j,k))
      ENDDO
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)*nb(3)
  ma_abs_r31 = 0
!----------------------
END FUNCTION ma_abs_r31
!===
INTEGER FUNCTION ma_chs_r31(nb,x,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(3),nbo,i,j,k,ij
  REAL :: x(nb(1),nb(2),nb(3)),y(nbo)
!---------------------------------------------------------------------
  ij = 0
  DO k=1,nb(3)
    DO j=1,nb(2)
      DO i=1,nb(1)
        ij = ij+1
        y(ij) = x(i,j,k)*(-1.)
      ENDDO
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)*nb(3)
  ma_chs_r31 = 0
!----------------------
END FUNCTION ma_chs_r31
!===
INTEGER FUNCTION ma_cels_r31(nb,x,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(3),nbo,i,j,k,ij
  REAL :: x(nb(1),nb(2),nb(3)),y(nbo)
!---------------------------------------------------------------------
  ij = 0
  DO k=1,nb(3)
    DO j=1,nb(2)
      DO i=1,nb(1)
        ij = ij+1
        y(ij) = x(i,j,k)-273.15
      ENDDO
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)*nb(3)
  ma_cels_r31 = 0
!-----------------------
END FUNCTION ma_cels_r31
!===
INTEGER FUNCTION ma_kelv_r31(nb,x,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(3),nbo,i,j,k,ij
  REAL :: x(nb(1),nb(2),nb(3)),y(nbo)
!---------------------------------------------------------------------
  ij = 0
  DO k=1,nb(3)
    DO j=1,nb(2)
      DO i=1,nb(1)
        ij = ij+1
        y(ij) = x(i,j,k)+273.15
      ENDDO
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)*nb(3)
  ma_kelv_r31 = 0
!-----------------------
  END FUNCTION ma_kelv_r31
!===
INTEGER FUNCTION ma_deg_r31(nb,x,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(3),nbo,i,j,k,ij
  REAL :: x(nb(1),nb(2),nb(3)),y(nbo)
!---------------------------------------------------------------------
  ij = 0
  DO k=1,nb(3)
    DO j=1,nb(2)
      DO i=1,nb(1)
        ij = ij+1
        y(ij) = x(i,j,k)*57.29577951
      ENDDO
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)*nb(3)
  ma_deg_r31 = 0
!----------------------
END FUNCTION ma_deg_r31
!===
INTEGER FUNCTION ma_rad_r31(nb,x,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(3),nbo,i,j,k,ij
  REAL :: x(nb(1),nb(2),nb(3)),y(nbo)
!---------------------------------------------------------------------
  ij = 0
  DO k=1,nb(3)
    DO j=1,nb(2)
      DO i=1,nb(1)
        ij = ij+1
        y(ij) = x(i,j,k)*0.01745329252
      ENDDO
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)*nb(3)
  ma_rad_r31 = 0
!----------------------
END FUNCTION ma_rad_r31
!===
INTEGER FUNCTION ma_ident_r31(nb,x,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(3),nbo,i,j,k,ij
  REAL :: x(nb(1),nb(2),nb(3)),y(nbo)
!---------------------------------------------------------------------
  ij = 0
  DO k=1,nb(3)
    DO j=1,nb(2)
      DO i=1,nb(1)
        ij = ij+1
        y(ij) = x(i,j,k)
      ENDDO
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)*nb(3)
  ma_ident_r31 = 0
!------------------------
END FUNCTION ma_ident_r31
!-
!=== OPERATIONS (two argument)
!-
INTEGER FUNCTION ma_add_r31(nb,x,s,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(3),nbo
  REAL :: x(nb(1),nb(2),nb(3)),s,y(nbo)
!-
  INTEGER :: i,j,k,ij
!---------------------------------------------------------------------
  ij = 0
  DO k=1,nb(3)
    DO j=1,nb(2)
      DO i=1,nb(1)
        ij = ij+1
        y(ij) = x(i,j,k)+s
      ENDDO
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)*nb(3)
  ma_add_r31 = 0
!----------------------
END FUNCTION ma_add_r31
!===
INTEGER FUNCTION ma_sub_r31(nb,x,s,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(3),nbo
  REAL :: x(nb(1),nb(2),nb(3)),s,y(nbo)
!-
  INTEGER :: i,j,k,ij
!---------------------------------------------------------------------
  ij = 0
  DO k=1,nb(3)
    DO j=1,nb(2)
      DO i=1,nb(1)
        ij = ij+1
        y(ij) = x(i,j,k)-s
      ENDDO
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)*nb(3)
  ma_sub_r31 = 0
!----------------------
END FUNCTION ma_sub_r31
!===
INTEGER FUNCTION ma_subi_r31(nb,x,s,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(3),nbo
  REAL :: x(nb(1),nb(2),nb(3)),s,y(nbo)
!-
  INTEGER :: i,j,k,ij
!---------------------------------------------------------------------
  ij = 0
  DO k=1,nb(3)
    DO j=1,nb(2)
      DO i=1,nb(1)
        ij = ij+1
        y(ij) =  s-x(i,j,k)
      ENDDO
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)*nb(3)
  ma_subi_r31 = 0
!-----------------------
END FUNCTION ma_subi_r31
!===
INTEGER FUNCTION ma_mult_r31(nb,x,s,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(3),nbo
  REAL :: x(nb(1),nb(2),nb(3)),s,y(nbo)
!-
  INTEGER :: i,j,k,ij
!---------------------------------------------------------------------
  ij = 0
  DO k=1,nb(3)
    DO j=1,nb(2)
      DO i=1,nb(1)
        ij = ij+1
        y(ij) = x(i,j,k)*s
      ENDDO
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)*nb(3)
  ma_mult_r31 = 0
!-----------------------
END FUNCTION ma_mult_r31
!===
INTEGER FUNCTION ma_div_r31(nb,x,s,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(3),nbo
  REAL :: x(nb(1),nb(2),nb(3)),s,y(nbo)
!-
  INTEGER :: i,j,k,ij
!---------------------------------------------------------------------
  ij = 0
  DO k=1,nb(3)
    DO j=1,nb(2)
      DO i=1,nb(1)
        ij = ij+1
        y(ij) = x(i,j,k)/s
      ENDDO
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)*nb(3)
  ma_div_r31 = 0
!----------------------
END FUNCTION ma_div_r31
!===
INTEGER FUNCTION ma_divi_r31(nb,x,s,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(3),nbo
  REAL :: x(nb(1),nb(2),nb(3)),s,y(nbo)
!-
  INTEGER :: i,j,k,ij
!---------------------------------------------------------------------
  ij = 0
  DO k=1,nb(3)
    DO j=1,nb(2)
      DO i=1,nb(1)
        ij = ij+1
        y(ij) = s/x(i,j,k)
      ENDDO
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)*nb(3)
  ma_divi_r31 = 0
!-----------------------
END FUNCTION ma_divi_r31
!===
INTEGER FUNCTION ma_power_r31(nb,x,s,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(3),nbo
  REAL :: x(nb(1),nb(2),nb(3)),s,y(nbo)
!-
  INTEGER :: i,j,k,ij
!---------------------------------------------------------------------
  ij = 0
  DO k=1,nb(3)
    DO j=1,nb(2)
      DO i=1,nb(1)
        ij = ij+1
        y(ij) = x(i,j,k)**s
      ENDDO
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)*nb(3)
  ma_power_r31 = 0
!------------------------
END FUNCTION ma_power_r31
!===
INTEGER FUNCTION ma_fumin_r31(nb,x,s,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(3),nbo
  REAL :: x(nb(1),nb(2),nb(3)),s,y(nbo)
!-
  INTEGER :: i,j,k,ij
!---------------------------------------------------------------------
  ij = 0
  DO k=1,nb(3)
    DO j=1,nb(2)
      DO i=1,nb(1)
        ij = ij+1
        y(ij) = MIN(x(i,j,k),s)
      ENDDO
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)*nb(3)
  ma_fumin_r31 = 0
!------------------------
END FUNCTION ma_fumin_r31
!===
INTEGER FUNCTION ma_fumax_r31(nb,x,s,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(3),nbo
  REAL :: x(nb(1),nb(2),nb(3)),s,y(nbo)
!-
  INTEGER :: i,j,k,ij
!---------------------------------------------------------------------
  ij = 0
  DO k=1,nb(3)
    DO j=1,nb(2)
      DO i=1,nb(1)
        ij = ij+1
        y(ij) = MAX(x(i,j,k),s)
      ENDDO
    ENDDO
  ENDDO
!-
  nbo = nb(1)*nb(2)*nb(3)
  ma_fumax_r31 = 0
!------------------------
END FUNCTION ma_fumax_r31
!===
INTEGER FUNCTION ma_fuscat_r31(nb,x,nbi,ind,miss_val,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(3),nbo,nbi
  INTEGER :: ind(nbi)
  REAL :: x(nb(1),nb(2),nb(3)),miss_val,y(nbo)
!-
  INTEGER :: i,j,k,ij,ii,ipos,ipp,isb
!---------------------------------------------------------------------
  ma_fuscat_r31 = 0
!-
  y(1:nbo) = miss_val
!-
  IF (nbi <= nb(1)*nb(2)*nb(3)) THEN
    ipos = 0
    isb = nb(1)*nb(2)
    DO ij=1,nbi
      IF (ind(ij) <= nbo .AND. ind(ij) > 0) THEN
        ipos = ipos+1
        k = ((ipos-1)/isb)+1
        ipp = ipos-(k-1)*isb
        j = ((ipp-1)/nb(1))+1
        i = (ipp-(j-1)*nb(1))
        y(ind(ij)) = x(i,j,k)
      ELSE
        IF (ind(ij) > nbo) ma_fuscat_r31  = ma_fuscat_r31+1
      ENDIF
    ENDDO
!-- Repeat the data if needed
    IF (MINVAL(ind) < 0) THEN
      DO i=1,nbi
        IF (ind(i) <= 0) THEN
          DO ii=1,ABS(ind(i))-1
            IF (ind(i+1)+ii <= nbo) THEN
              y(ind(i+1)+ii) = y(ind(i+1))
            ELSE
              ma_fuscat_r31  = ma_fuscat_r31+1
            ENDIF
          ENDDO
        ENDIF
      ENDDO
    ENDIF
  ELSE
    ma_fuscat_r31  = 1
  ENDIF
!-------------------------
END FUNCTION ma_fuscat_r31
!===
INTEGER FUNCTION ma_fugath_r31(nb,x,nbi,ind,miss_val,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(3),nbo,nbi
  INTEGER :: ind(nbi)
  REAL :: x(nb(1),nb(2),nb(3)),miss_val,y(nbo)
!-
  INTEGER :: i,j,k,ij,ipos,ipp,isb
!---------------------------------------------------------------------
  IF (nbi <= nbo) THEN
    ma_fugath_r31 = 0
    y(1:nbo) = miss_val
    ipos = 0
    isb = nb(1)*nb(2)
    DO ij=1,nbi
      IF (ipos+1 <= nbo) THEN
        IF (ind(ij) > 0) THEN
          k = ((ind(ij)-1)/isb)+1
          ipp = ind(ij)-(k-1)*isb
          j = ((ipp-1)/nb(1))+1
          i = (ipp-(j-1)*nb(1))
          ipos = ipos+1
          y(ipos) = x(i,j,k)
        ENDIF
      ELSE
        IF (ipos+1 > nbo) ma_fugath_r31  = ma_fugath_r31+1
      ENDIF
    ENDDO
  ELSE
    ma_fugath_r31 = 1
  ENDIF
  nbo = ipos
!-------------------------
END FUNCTION ma_fugath_r31
!===
INTEGER FUNCTION ma_fufill_r31(nb,x,nbi,ind,miss_val,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(3),nbo,nbi
  INTEGER :: ind(nbi)
  REAL :: x(nb(1),nb(2),nb(3)),miss_val,y(nbo)
!-
  INTEGER :: i,j,k,ij,ii,ipos,ipp,isb
!---------------------------------------------------------------------
  ma_fufill_r31 = 0
  IF (nbi <= nb(1)*nb(2)*nb(3)) THEN
    ipos = 0
    isb = nb(1)*nb(2)
    DO ij=1,nbi
      IF (ind(ij) <= nbo .AND. ind(ij) > 0) THEN
        ipos = ipos+1
        k = ((ipos-1)/isb)+1
        ipp = ipos-(k-1)*isb
        j = ((ipp-1)/nb(1))+1
        i = (ipp-(j-1)*nb(1))
        y(ind(ij)) = x(i,j,k)
      ELSE
        IF (ind(ij) > nbo) ma_fufill_r31  = ma_fufill_r31+1
      ENDIF
    ENDDO
!-- Repeat the data if needed
    IF (MINVAL(ind) < 0) THEN
      DO i=1,nbi
        IF (ind(i) <= 0) THEN
          DO ii=1,ABS(ind(i))-1
            IF (ind(i+1)+ii <= nbo) THEN
              y(ind(i+1)+ii) = y(ind(i+1))
            ELSE
              ma_fufill_r31  = ma_fufill_r31+1
            ENDIF
          ENDDO
        ENDIF
      ENDDO
    ENDIF
  ELSE
    ma_fufill_r31  = 1
  ENDIF
!-------------------------
END FUNCTION ma_fufill_r31
!===
INTEGER FUNCTION ma_fucoll_r31(nb,x,nbi,ind,miss_val,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(3),nbo,nbi
  INTEGER :: ind(nbi)
  REAL :: x(nb(1),nb(2),nb(3)),miss_val,y(nbo)
!-
  INTEGER :: i,j,k,ij,ipos,ipp,isb
!---------------------------------------------------------------------
  IF (nbi <= nbo) THEN
    ma_fucoll_r31 = 0
    ipos = 0
    isb = nb(1)*nb(2)
    DO ij=1,nbi
      IF (ipos+1 <= nbo) THEN
        IF (ind(ij) > 0) THEN
          k = ((ind(ij)-1)/isb)+1
          ipp = ind(ij)-(k-1)*isb
          j = ((ipp-1)/nb(1))+1
          i = (ipp-(j-1)*nb(1))
          ipos = ipos+1
          y(ipos) = x(i,j,k)
        ENDIF
      ELSE
        IF (ipos+1 > nbo) ma_fucoll_r31  = ma_fucoll_r31+1
      ENDIF
    ENDDO
  ELSE
    ma_fucoll_r31 = 1
  ENDIF
  nbo = ipos
!-------------------------
END FUNCTION ma_fucoll_r31
!===
INTEGER FUNCTION ma_fuundef_r31(nb,x,nbi,ind,miss_val,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(3),nbo,nbi
  INTEGER :: ind(nbi)
  REAL :: x(nb(1),nb(2),nb(3)),miss_val,y(nbo)
!-
  INTEGER :: i,j,k,ij,ipp,isb
!---------------------------------------------------------------------
  IF (nbi <= nbo .AND. nbo == nb(1)*nb(2)*nb(3)) THEN
    ma_fuundef_r31 = 0
    isb = nb(1)*nb(2)
    DO ij=1,nbo
      k = ((ij-1)/isb)+1
      ipp = ij-(k-1)*isb
      j = ((ipp-1)/nb(1))+1
      i = (ipp-(j-1)*nb(1))
      y(ij) = x(i,j,k)
    ENDDO
    DO i=1,nbi
      IF (ind(i) <= nbo .AND. ind(i) > 0) THEN
        y(ind(i)) =  miss_val
      ELSE
        IF (ind(i) > nbo) ma_fuundef_r31  = ma_fuundef_r31+1
      ENDIF
    ENDDO
  ELSE
    ma_fuundef_r31 = 1
  ENDIF
!--------------------------
END FUNCTION ma_fuundef_r31
!===
INTEGER FUNCTION ma_fuonly_r31(nb,x,nbi,ind,miss_val,nbo,y)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  INTEGER :: nb(3),nbo,nbi
  INTEGER :: ind(nbi)
  REAL :: x(nb(1),nb(2),nb(3)),miss_val,y(nbo)
!-
  INTEGER :: i,j,k,ij,ipp,isb
!---------------------------------------------------------------------
  IF (     (nbi <= nbo).AND.(nbo == nb(1)*nb(2)*nb(3)) &
 &    .AND.ALL(ind(1:nbi) <= nbo) ) THEN
    ma_fuonly_r31 = 0
    y(1:nbo) = miss_val
    isb = nb(1)*nb(2)
    DO ij=1,nbi
      IF (ind(ij) > 0) THEN
        k = ((ind(ij)-1)/isb)+1
        ipp = ind(ij)-(k-1)*isb
        j = ((ipp-1)/nb(1))+1
        i = (ipp-(j-1)*nb(1))
        y(ind(ij)) = x(i,j,k)
      ENDIF
    ENDDO
  ELSE
    ma_fuonly_r31 = 1
  ENDIF
!-------------------------
END FUNCTION ma_fuonly_r31
!===
SUBROUTINE moycum (opp,np,px,py,pwx)
!---------------------------------------------------------------------
!- Does time operations
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  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 ipslerr(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 moycum_index( opp, px, py, pwx, nbi, ind )
!---------------------------------------------------------------------
!- Does time operations on index points 
!---------------------------------------------------------------------
  IMPLICIT NONE
!-

  !! 0. Parameters and variables declaration

  !! 0.1 Input variables

  CHARACTER(LEN=7), INTENT(in)        :: opp  !! Operation performed
  INTEGER, INTENT(in)                 :: nbi  !! Size of index vector
  INTEGER, DIMENSION(nbi), INTENT(in) :: ind  !! Index vector
  REAL, DIMENSION(:), INTENT(in)      :: py   !! Vector containing the 
                                              !! previous values of px
                                              !! Warning : due to memory 
                                              !! optimization, we have
                                              !! generally SIZE(px) /= SIZE(py)  
  INTEGER, INTENT(in)                 :: pwx  !! Used to calculate average value                             

  !! 0.3 Modified variables

  REAL, DIMENSION(:), INTENT(inout)   :: px   !! Result

  !! 0.4 Local variables

  INTEGER :: ig                               !! Index

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

  !! Perform operations only if the values of ind don't exceed the size of px 

  IF ( MAXVAL(ind) > SIZE(px) ) THEN
     CALL ipslerr(3,"moycum_index", &
          & "the index vector is out of range for px", &
          & "Indexation vector problem. We stop", " " )
  END IF

  IF (pwx /= 0) THEN
     IF      (opp == 'ave') THEN
        DO ig = 1,nbi
           px(ind(ig)) = (px(ind(ig))*pwx + py(ind(ig)))/REAL(pwx+1)
        END DO
     ELSE IF (opp == 't_sum') THEN
        DO ig = 1,nbi
           px(ind(ig)) = px(ind(ig)) + py(ind(ig))
        END DO
     ELSE IF ( (opp == 'l_min').OR.(opp == 't_min') ) THEN
        DO ig = 1,nbi
           px(ind(ig)) = MIN(px(ind(ig)),py(ind(ig)))
        END DO
     ELSE IF ( (opp == 'l_max').OR.(opp == 't_max') ) THEN
        DO ig = 1,nbi
           px(ind(ig)) = MAX(px(ind(ig)),py(ind(ig)))
        END DO
     ELSE
        CALL ipslerr(3,"moycum_index",'Unknown time operation',opp,' ')
     END IF
  ELSE
    IF      (opp == 'l_min') THEN
       DO ig = 1,nbi
          px(ind(ig)) = MIN(px(ind(ig)),py(ind(ig)))
       END DO
    ELSE IF (opp == 'l_max') THEN
       DO ig = 1,nbi
          px(ind(ig)) = MAX(px(ind(ig)),py(ind(ig)))
       END DO
    ELSE
       DO ig = 1,nbi
          px(ind(ig)) = py(ind(ig))
       END DO
    ENDIF
 END IF

END SUBROUTINE moycum_index


!-----------------
END MODULE mathelp
