source: IOIPSL/trunk/src/mathelp.f90 @ 3279

Last change on this file since 3279 was 1927, checked in by dsolyga, 11 years ago

Introduced the new subroutine moycum_index. Works the same way as moycum but make computations only on index points. Used only when scatter operation is performed. Help to reduce the computational time of ORCHIDEE. For the other models, should not change the results.

  • Property svn:keywords set to Id
File size: 82.9 KB
Line 
1MODULE mathelp
2!-
3!$Id$
4!-
5! This software is governed by the CeCILL license
6! See IOIPSL/IOIPSL_License_CeCILL.txt
7!---------------------------------------------------------------------
8  USE errioipsl,ONLY : ipslerr
9  USE stringop
10!-
11  PRIVATE
12  PUBLIC :: mathop,moycum,moycum_index,buildop
13!-
14  INTERFACE mathop
15    MODULE PROCEDURE mathop_r11,mathop_r21,mathop_r31
16  END INTERFACE
17!-
18!- Variables used to detect and identify the operations
19!-
20  CHARACTER(LEN=80),SAVE :: &
21 &  seps='( ) , + - / * ^', ops = '+ - * / ^', mima = 'min max'
22  CHARACTER(LEN=250),SAVE :: &
23 &  funcs = 'sin cos tan asin acos atan exp log sqrt chs abs '&
24 & //'cels kelv deg rad gather scatter fill coll undef only ident'
25  CHARACTER(LEN=120),SAVE :: &
26 &  indexfu = 'gather, scatter, fill, coll, undef, only'
27!---------------------------------------------------------------------
28CONTAINS
29!===
30SUBROUTINE buildop (c_str,ex_topps,topp,fill_val,opps,scal,nbops)
31!---------------------------------------------------------------------
32!- This subroutine decomposes the input string in the elementary
33!- functions which need to be applied to the vector of data.
34!- This vector is represented by X in the string.
35!- This subroutine is the driver of the decomposition and gets
36!- the time operation but then call decoop for the other operations
37!-
38!- INPUT
39!-
40!- c_str    : String containing the operations
41!- ex_toops : Time operations that can be expected within the string
42!- fill_val :
43!-
44!- OUTPUT
45!-
46!- topp   : Time operation
47!- opps   :
48!- scal   :
49!- nbops  :
50!---------------------------------------------------------------------
51  IMPLICIT NONE
52!-
53  CHARACTER(LEN=*),INTENT(IN) :: c_str,ex_topps
54  CHARACTER(LEN=*),INTENT(OUT) :: topp
55  CHARACTER(LEN=*),DIMENSION(:),INTENT(OUT) :: opps
56  REAL,INTENT(IN) :: fill_val
57  REAL,DIMENSION(:),INTENT(OUT) :: scal
58  INTEGER,INTENT(OUT) :: nbops
59!-
60  CHARACTER(LEN=LEN(c_str)) :: str,new_str
61  INTEGER :: leng,ind_opb,ind_clb
62!-
63  LOGICAL :: check = .FALSE.
64!---------------------------------------------------------------------
65  IF (check) WRITE(*,*) 'buildop : Some preliminary cleaning'
66!-
67  str = c_str
68  leng = LEN_TRIM(str)
69  IF ( str(1:1) == '(' .AND. str(leng:leng) == ')' ) THEN
70    str = str(2:leng-1)
71    leng = leng-2
72  ENDIF
73!-
74  IF (check) &
75 &  WRITE(*,*) 'buildop : Starting to test the various options'
76!-
77  IF (leng <= 5 .AND. INDEX(ex_topps,str(1:leng)) > 0) THEN
78    IF (check) WRITE(*,*) 'buildop : Time operation only'
79    nbops = 0
80    topp = str(1:leng)
81  ELSE
82    IF (check) THEN
83      WRITE(*,*) 'buildop : Time operation and something else'
84    ENDIF
85!--
86    ind_opb = INDEX(str(1:leng),'(')
87    IF (ind_opb > 0) THEN
88      IF (INDEX(ex_topps,str(1:ind_opb-1)) > 0) THEN
89        IF (check) THEN
90          WRITE(*,'(2a)') &
91 &          ' buildop : Extract time operation from : ',str
92        ENDIF
93        topp = str(1:ind_opb-1)
94        ind_clb = INDEX(str(1:leng),')',BACK=.TRUE.)
95        new_str = str(ind_opb+1:ind_clb-1)
96        IF (check) THEN
97          WRITE(*,'(2a,2I3)') &
98 &          ' buildop : Call decoop ',new_str,ind_opb,ind_clb
99        ENDIF
100        CALL decoop (new_str,fill_val,opps,scal,nbops)
101      ELSE
102        CALL ipslerr(3,'buildop', &
103 &        'time operation does not exist',str(1:ind_opb-1),' ')
104      ENDIF
105    ELSE
106      CALL ipslerr(3,'buildop', &
107 &      'some long operation exists but wihout parenthesis', &
108 &      str(1:leng),' ')
109    ENDIF
110  ENDIF
111!-
112  IF (check) THEN
113    DO leng=1,nbops
114      WRITE(*,*) &
115 &      'buildop : i -- opps, scal : ',leng,opps(leng),scal(leng)
116    ENDDO
117  ENDIF
118!---------------------
119END SUBROUTINE buildop
120!===
121SUBROUTINE decoop (pstr,fill_val,opps,scal,nbops)
122!---------------------------------------------------------------------
123  IMPLICIT NONE
124!-
125  CHARACTER(LEN=*),INTENT(IN) :: pstr
126  REAL,INTENT(IN) :: fill_val
127  CHARACTER(LEN=*),DIMENSION(:),INTENT(OUT) :: opps
128  REAL,DIMENSION(:),INTENT(OUT) :: scal
129  INTEGER,INTENT(OUT) :: nbops
130!-
131  CHARACTER(LEN=1),DIMENSION(2) :: f_char,s_char
132  INTEGER,DIMENSION(2) :: f_pos,s_pos
133  CHARACTER(LEN=20) :: opp_str,scal_str
134  CHARACTER(LEN=LEN(pstr)) :: str
135  INTEGER :: nbsep,nbops_max,xpos,leng,ppos,epos,int_tmp
136  CHARACTER(LEN=3) :: tl,dl
137  CHARACTER(LEN=10) :: fmt
138!-
139  LOGICAL :: check = .FALSE.,prio
140!---------------------------------------------------------------------
141  IF (check) WRITE(*,'(2A)') ' decoop : Incoming string : ',pstr
142!-
143  str = pstr; nbops = 0;
144!-
145  CALL findsep (str,nbsep,f_char,f_pos,s_char,s_pos)
146  IF (check) WRITE(*,*) 'decoop : Out of findsep',nbsep
147!-
148  nbops_max = min(SIZE(opps),SIZE(scal))
149!-
150  DO WHILE (nbsep > 0)
151    IF (nbops >= nbops_max) THEN
152      CALL ipslerr(3,'decoop','Expression too complex',TRIM(str),' ')
153    ENDIF
154!--
155    xpos = INDEX(str,'X')
156    leng = LEN_TRIM(str)
157    nbops = nbops+1
158!--
159    IF (check) THEN
160      WRITE(*,*) 'decoop : str   -> ',TRIM(str)
161      WRITE(*,*) 'decoop : nbops -> ',nbops
162      WRITE(*,*) s_char(1),'-',f_char(1),'|',f_char(2),'-',s_char(2)
163      WRITE(*,*) s_pos(1),'-',f_pos(1),'|',f_pos(2),'-',s_pos(2)
164    ENDIF
165!---
166!-- Start the analysis of the syntax. 3 types of constructs
167!-- are recognized.  They are scanned sequentialy
168!---
169    IF (nbsep == 1) THEN
170      IF (check) WRITE(*,*) 'decoop : Only one operation'
171      IF (INDEX(ops,f_char(1)) > 0) THEN
172!------ Type : scal+X
173        IF (f_char(1) == '-' .OR. f_char(1) == '/') THEN
174          opp_str = f_char(1)//'I'
175        ELSE
176          opp_str = f_char(1)
177        ENDIF
178        scal_str = str(s_pos(1)+1:f_pos(1)-1)
179        str = 'X'
180      ELSE IF (INDEX(ops,f_char(2)) > 0) THEN
181!------ Type : X+scal
182        opp_str = f_char(2)
183        scal_str = str(f_pos(2)+1:s_pos(2)-1)
184        str = 'X'
185      ELSE
186        CALL ipslerr(3,'decoop', &
187 &        'Unknown operations of type X+scal',f_char(1),pstr)
188      ENDIF
189    ELSE
190      IF (check) WRITE(*,*) 'decoop : More complex operation'
191      IF ( f_char(1) == '(' .AND. f_char(2) == ')' ) THEN
192!------ Type : sin(X)
193        opp_str = str(s_pos(1)+1:f_pos(1)-1)
194        scal_str = '?'
195        str = str(1:s_pos(1))//'X'//str(f_pos(2)+1:leng)
196      ELSE IF (    (f_char(1) == '(' .AND. f_char(2) == ',')&
197 &             .OR.(f_char(1) == ',' .AND. f_char(2) == ')')) THEN
198!------ Type : max(X,scal) or max(scal,X)
199        IF (f_char(1) == '(' .AND. s_char(2) == ')') THEN
200!-------- Type : max(X,scal)
201          opp_str = str(f_pos(1)-3:f_pos(1)-1)
202          scal_str = str(f_pos(2)+1:s_pos(2)-1)
203          str = str(1:f_pos(1)-4)//'X'//str(s_pos(2)+1:leng)
204        ELSE IF (f_char(1) == ',' .AND. s_char(1) == '(') THEN
205!-------- Type : max(scal,X)
206          opp_str = str(s_pos(1)-3:s_pos(1)-1)
207          scal_str = str(s_pos(1)+1:f_pos(1)-1)
208          str = str(1:s_pos(1)-4)//'X'//str(f_pos(2)+1:leng)
209        ELSE
210          CALL ipslerr(3,'decoop','Syntax error 1',str,' ')
211        ENDIF
212      ELSE
213        prio = (f_char(2) == '*').OR.(f_char(2) == '^')
214        IF (     (INDEX(ops,f_char(1)) > 0) &
215 &          .AND.(xpos-f_pos(1) == 1).AND.(.NOT.prio) ) THEN
216!-------- Type : ... scal+X ...
217          IF (f_char(1) == '-' .OR. f_char(1) == '/') THEN
218            opp_str = f_char(1)//'I'
219          ELSE
220            opp_str = f_char(1)
221          ENDIF
222          scal_str = str(s_pos(1)+1:f_pos(1)-1)
223          str = str(1:s_pos(1))//'X'//str(f_pos(1)+2:leng)
224        ELSE IF (     (INDEX(ops,f_char(2)) > 0) &
225 &               .AND.(f_pos(2)-xpos == 1) ) THEN
226!-------- Type : ... X+scal ...
227          opp_str = f_char(2)
228          scal_str = str(f_pos(2)+1:s_pos(2)-1)
229          str = str(1:f_pos(2)-2)//'X'//str(s_pos(2):leng)
230        ELSE
231          CALL ipslerr(3,'decoop','Syntax error 2',str,' ')
232        ENDIF
233      ENDIF
234    ENDIF
235!---
236    IF (check) WRITE(*,*) 'decoop : Finished syntax,str = ',TRIM(str)
237!---
238!-- Now that the different components of the operation are identified
239!-- we transform them into what is going to be used in the program
240!---
241    IF (INDEX(scal_str,'?') > 0) THEN
242      IF (INDEX(funcs,opp_str(1:LEN_TRIM(opp_str))) > 0) THEN
243        opps(nbops) = opp_str(1:LEN_TRIM(opp_str))
244        scal(nbops) = fill_val
245      ELSE
246        CALL ipslerr(3,'decoop', &
247 &        'Unknown function',opp_str(1:LEN_TRIM(opp_str)),' ')
248      ENDIF
249    ELSE
250      leng = LEN_TRIM(opp_str)
251      IF (INDEX(mima,opp_str(1:leng)) > 0) THEN
252        opps(nbops) = 'fu'//opp_str(1:leng)
253      ELSE
254        IF (INDEX(opp_str(1:leng),'+') > 0) THEN
255          opps(nbops) = 'add'
256        ELSE IF (INDEX(opp_str(1:leng),'-I') > 0) THEN
257          opps(nbops) = 'subi'
258        ELSE IF (INDEX(opp_str(1:leng),'-') > 0) THEN
259          opps(nbops) = 'sub'
260        ELSE IF (INDEX(opp_str(1:leng),'*') > 0) THEN
261          opps(nbops) = 'mult'
262        ELSE IF (INDEX(opp_str(1:leng),'/') > 0) THEN
263          opps(nbops) = 'div'
264        ELSE IF (INDEX(opp_str(1:leng),'/I') > 0) THEN
265          opps(nbops) = 'divi'
266        ELSE IF (INDEX(opp_str(1:leng),'^') > 0) THEN
267          opps(nbops) = 'power'
268        ELSE
269          CALL ipslerr(3,'decoop', &
270 &          'Unknown operation',opp_str(1:leng),' ')
271        ENDIF
272      ENDIF
273!-----
274      leng = LEN_TRIM(scal_str)
275      ppos = INDEX(scal_str,'.')
276      epos = INDEX(scal_str,'e')
277      IF (epos == 0) epos = INDEX(scal_str,'E')
278!-----
279!---- Try to catch a few errors
280!-----
281      IF (INDEX(ops,scal_str) > 0) THEN
282        CALL ipslerr(3,'decoop', &
283 &        'Strange scalar you have here ',scal_str,pstr)
284      ENDIF
285      IF (epos > 0) THEN
286        WRITE(tl,'(I3.3)') leng
287        WRITE(dl,'(I3.3)') epos-ppos-1
288        fmt='(e'//tl//'.'//dl//')'
289        READ(scal_str,fmt) scal(nbops)
290      ELSE IF (ppos > 0) THEN
291        WRITE(tl,'(I3.3)') leng
292        WRITE(dl,'(I3.3)') leng-ppos
293        fmt='(f'//tl//'.'//dl//')'
294        READ(scal_str,fmt) scal(nbops)
295      ELSE
296        WRITE(tl,'(I3.3)') leng
297        fmt = '(I'//tl//')'
298        READ(scal_str,fmt) int_tmp
299        scal(nbops) = REAL(int_tmp)
300      ENDIF
301    ENDIF
302    IF (check) WRITE(*,*) 'decoop : Finished interpretation'
303    CALL findsep(str,nbsep,f_char,f_pos,s_char,s_pos)
304  ENDDO
305!--------------------
306END SUBROUTINE decoop
307!===
308SUBROUTINE findsep (str,nbsep,f_char,f_pos,s_char,s_pos)
309!---------------------------------------------------------------------
310!- Subroutine finds all separators in a given string
311!- It returns the following information about str :
312!-   f_char : The first separation character
313!-            (1 for before and 2 for after)
314!-   f_pos  : The position of the first separator
315!-   s_char : The second separation character
316!-            (1 for before and 2 for after)
317!-   s_pos  : The position of the second separator
318!---------------------------------------------------------------------
319  IMPLICIT NONE
320!-
321  CHARACTER(LEN=*),INTENT(INOUT) :: str
322  INTEGER :: nbsep
323  CHARACTER(LEN=1),DIMENSION(2) :: f_char,s_char
324  INTEGER,DIMENSION(2) :: f_pos,s_pos
325!-
326  CHARACTER(LEN=10) :: str_tmp
327  LOGICAL :: f_found,s_found
328  INTEGER :: ind,xpos,leng,i
329!-
330  LOGICAL :: check = .FALSE.
331!---------------------------------------------------------------------
332  IF (check) WRITE(*,*) 'findsep : call cleanstr: ',TRIM(str)
333!-
334  CALL cleanstr(str)
335!-
336  IF (check) WRITE(*,*) 'findsep : out of cleanstr: ',TRIM(str)
337!-
338  xpos = INDEX(str,'X')
339  leng = LEN_TRIM(str)
340!-
341  f_pos(1:2) = (/ 0,leng+1 /)
342  f_char(1:2) = (/ '?','?' /)
343  s_pos(1:2) = (/ 0,leng+1 /)
344  s_char(1:2) = (/ '?','?' /)
345!-
346  nbsep = 0
347!-
348  f_found = .FALSE.
349  s_found = .FALSE.
350  IF (xpos > 1) THEN
351    DO i=xpos-1,1,-1
352      ind = INDEX(seps,str(i:i))
353      IF (ind > 0) THEN
354        IF (.NOT.f_found) THEN
355          f_char(1) = str(i:i)
356          f_pos(1) = i
357          nbsep = nbsep+1
358          f_found = .TRUE.
359        ELSE IF (.NOT.s_found) THEN
360          s_char(1) = str(i:i)
361          s_pos(1) = i
362          nbsep = nbsep+1
363          s_found = .TRUE.
364        ENDIF
365      ENDIF
366    ENDDO
367  ENDIF
368!-
369  f_found = .FALSE.
370  s_found = .FALSE.
371  IF (xpos < leng) THEN
372    DO i=xpos+1,leng
373      ind = INDEX(seps,str(i:i))
374      IF (ind > 0) THEN
375        IF (.NOT.f_found) THEN
376          f_char(2) = str(i:i)
377          f_pos(2) = i
378          nbsep = nbsep+1
379          f_found = .TRUE.
380        ELSE IF (.NOT.s_found) THEN
381          s_char(2) = str(i:i)
382          s_pos(2) = i
383          nbsep = nbsep+1
384          s_found = .TRUE.
385        ENDIF
386      ENDIF
387    ENDDO
388  ENDIF
389!-
390  IF (nbsep > 4) THEN
391    WRITE(str_tmp,'("number :",I3)') nbsep
392    CALL ipslerr(3,'findsep', &
393 &    'How can I find that many separators',str_tmp,TRIM(str))
394  ENDIF
395!-
396  IF (check) WRITE(*,*) 'Finished findsep : ',nbsep,leng
397!---------------------
398END SUBROUTINE findsep
399!===
400SUBROUTINE cleanstr(str)
401!---------------------------------------------------------------------
402!- We clean up the string by taking out the extra () and puting
403!- everything in lower case except for the X describing the variable
404!---------------------------------------------------------------------
405  IMPLICIT NONE
406!-
407  CHARACTER(LEN=*),INTENT(INOUT) :: str
408!-
409  INTEGER :: ind,leng,ic,it
410  LOGICAL :: check = .FALSE.
411!---------------------------------------------------------------------
412  leng = LEN_TRIM(str)
413  CALL strlowercase(str)
414!-
415  ind = INDEX(str,'x')
416  IF (check) THEN
417    WRITE (*,*) 'cleanstr 1.0 : ind = ',ind, &
418&               ' str = ',str(1:leng),'---'
419  ENDIF
420!-
421! If the character before the x is not a letter then we can assume
422! that it is the variable and promote it to a capital letter
423!-
424  DO WHILE (ind > 0)
425    ic = 0
426    IF (ind > 1) ic = IACHAR(str(ind-1:ind-1))
427    IF (ic < 97 .OR. ic > 122) THEN
428      str(ind:ind) = 'X'
429    ENDIF
430    it = INDEX(str(ind+1:leng),'x')
431    IF (it > 0) THEN
432      ind = ind+it
433    ELSE
434      ind = it
435    ENDIF
436  ENDDO
437!-
438  IF (check) WRITE (*,*) 'cleanstr 2.0 : str = ',str(1:leng),'---'
439!-
440  IF ( str(1:1) == '(' .AND. str(leng:leng) == ')' ) THEN
441    str = str(2:leng-1)
442  ENDIF
443!-
444  IF (check) WRITE (*,*) 'cleanstr 3.0 : str = ',str(1:leng),'---'
445!-
446  leng = LEN_TRIM(str)
447  ind = INDEX(str,'((X))')
448  IF (ind > 0) THEN
449    str=str(1:ind-1)//'(X)'//str(ind+5:leng)//'  '
450  ENDIF
451!-
452  IF (check) WRITE (*,*) 'cleanstr 4.0 : str = ',str(1:leng),'---'
453!-
454  leng = LEN_TRIM(str)
455  ind = INDEX(str,'(X)')
456  IF (ind > 0 .AND. ind+3 < leng) THEN
457    IF (      (INDEX(seps,str(ind-1:ind-1)) > 0) &
458 &      .AND. (INDEX(seps,str(ind+3:ind+3)) > 0) ) THEN
459      str=str(1:ind-1)//'X'//str(ind+3:leng)//'  '
460    ENDIF
461  ENDIF
462!-
463  IF (check) WRITE (*,*) 'cleanstr 5.0 : str = ',str(1:leng),'---'
464!-
465  leng = LEN_TRIM(str)
466  ind = INDEX(str(1:leng),' ')
467  DO WHILE (ind > 0)
468    str=str(1:ind-1)//str(ind+1:leng)//' '
469    leng = LEN_TRIM(str)
470    ind = INDEX(str(1:leng),' ')
471  ENDDO
472!-
473  IF (check) WRITE (*,*) 'cleanstr 6.0 : str = ',str(1:leng),'---'
474!----------------------
475END SUBROUTINE cleanstr
476!===
477!===
478SUBROUTINE mathop_r11 &
479 &  (fun,nb,work_in,miss_val,nb_index,nindex,scal,nb_max,work_out)
480!---------------------------------------------------------------------
481!- This subroutines gives an interface to the various operation
482!- which are allowed. The interface is general enough to allow its use
483!- for other cases.
484!-
485!- INPUT
486!-
487!- fun      : function to be applied to the vector of data
488!- nb       : Length of input vector
489!- work_in  : Input vector of data (REAL)
490!- miss_val : The value of the missing data flag (it has to be a
491!-            maximum value, in f90 : huge( a real ))
492!- nb_index : Length of index vector
493!- nindex   : Vector of indices
494!- scal     : A scalar value for vector/scalar operations
495!- nb_max   : maximum length of output vector
496!-
497!- OUTPUT
498!-
499!- nb_max   : Actual length of output variable
500!- work_out : Output vector after the operation was applied
501!---------------------------------------------------------------------
502  IMPLICIT NONE
503!-
504  CHARACTER(LEN=7) :: fun
505  INTEGER :: nb,nb_max,nb_index
506  INTEGER :: nindex(nb_index)
507  REAL :: work_in(nb),scal,miss_val
508  REAL :: work_out(nb_max)
509!-
510  INTEGER :: ierr
511!-
512  INTRINSIC SIN,COS,TAN,ASIN,ACOS,ATAN,EXP,LOG,SQRT,ABS
513!---------------------------------------------------------------------
514  ierr = 0
515!-
516  IF (scal >= miss_val-1.) THEN
517    IF (INDEX(indexfu,fun(1:LEN_TRIM(fun))) == 0) THEN
518      SELECT CASE (fun)
519      CASE('sin')
520        ierr = ma_sin_r11(nb,work_in,nb_max,work_out)
521      CASE('cos')
522        ierr = ma_cos_r11(nb,work_in,nb_max,work_out)
523      CASE('tan')
524        ierr = ma_tan_r11(nb,work_in,nb_max,work_out)
525      CASE('asin')
526        ierr = ma_asin_r11(nb,work_in,nb_max,work_out)
527      CASE('acos')
528        ierr = ma_acos_r11(nb,work_in,nb_max,work_out)
529      CASE('atan')
530        ierr = ma_atan_r11(nb,work_in,nb_max,work_out)
531      CASE('exp')
532        ierr = ma_exp_r11(nb,work_in,nb_max,work_out)
533      CASE('log')
534        ierr = ma_log_r11(nb,work_in,nb_max,work_out)
535      CASE('sqrt')
536        ierr = ma_sqrt_r11(nb,work_in,nb_max,work_out)
537      CASE('chs')
538        ierr = ma_chs_r11(nb,work_in,nb_max,work_out)
539      CASE('abs')
540        ierr = ma_abs_r11(nb,work_in,nb_max,work_out)
541      CASE('cels')
542        ierr = ma_cels_r11(nb,work_in,nb_max,work_out)
543      CASE('kelv')
544        ierr = ma_kelv_r11(nb,work_in,nb_max,work_out)
545      CASE('deg')
546        ierr = ma_deg_r11(nb,work_in,nb_max,work_out)
547      CASE('rad')
548        ierr = ma_rad_r11(nb,work_in,nb_max,work_out)
549      CASE('ident')
550        ierr = ma_ident_r11(nb,work_in,nb_max,work_out)
551      CASE DEFAULT
552        CALL ipslerr(3,"mathop", &
553 &        'scalar variable undefined and no indexing', &
554 &        'but still unknown function',fun)
555      END SELECT
556      IF (ierr > 0) THEN
557        CALL ipslerr(3,"mathop", &
558 &        'Error while executing a simple function',fun,' ')
559      ENDIF
560    ELSE
561      SELECT CASE (fun)
562      CASE('gather')
563        ierr = ma_fugath_r11(nb,work_in,nb_index,nindex, &
564&                            miss_val,nb_max,work_out)
565      CASE('scatter')
566        IF (nb_index > nb) THEN
567          work_out(1:nb_max) = miss_val
568          ierr=1
569        ELSE
570          ierr = ma_fuscat_r11(nb,work_in,nb_index,nindex, &
571&                              miss_val,nb_max,work_out)
572        ENDIF
573      CASE('coll')
574        ierr = ma_fucoll_r11(nb,work_in,nb_index,nindex, &
575&                            miss_val,nb_max,work_out)
576      CASE('fill')
577        ierr = ma_fufill_r11(nb,work_in,nb_index,nindex, &
578&                            miss_val,nb_max,work_out)
579      CASE('undef')
580        ierr = ma_fuundef_r11(nb,work_in,nb_index,nindex, &
581&                            miss_val,nb_max,work_out)
582      CASE('only')
583        ierr = ma_fuonly_r11(nb,work_in,nb_index,nindex, &
584&                            miss_val,nb_max,work_out)
585      CASE DEFAULT
586        CALL ipslerr(3,"mathop", &
587 &        'scalar variable undefined and indexing',&
588 &        'was requested but with unknown function',fun)
589      END SELECT
590      IF (ierr > 0) THEN
591        CALL ipslerr(3,"mathop_r11", &
592 &        'Error while executing an indexing function',fun,' ')
593      ENDIF
594    ENDIF
595  ELSE
596    SELECT CASE (fun)
597    CASE('fumin')
598      ierr = ma_fumin_r11(nb,work_in,scal,nb_max,work_out)
599    CASE('fumax')
600      ierr = ma_fumax_r11(nb,work_in,scal,nb_max,work_out)
601    CASE('add')
602      ierr = ma_add_r11(nb,work_in,scal,nb_max,work_out)
603    CASE('subi')
604      ierr = ma_subi_r11(nb,work_in,scal,nb_max,work_out)
605    CASE('sub')
606      ierr = ma_sub_r11(nb,work_in,scal,nb_max,work_out)
607    CASE('mult')
608      ierr = ma_mult_r11(nb,work_in,scal,nb_max,work_out)
609    CASE('div')
610      ierr = ma_div_r11(nb,work_in,scal,nb_max,work_out)
611    CASE('divi')
612      ierr = ma_divi_r11(nb,work_in,scal,nb_max,work_out)
613    CASE('power')
614      ierr = ma_power_r11(nb,work_in,scal,nb_max,work_out)
615    CASE DEFAULT
616      CALL ipslerr(3,"mathop", &
617 &      'Unknown operation with a scalar',fun,' ')
618    END SELECT
619    IF (ierr > 0) THEN
620      CALL ipslerr(3,"mathop", &
621 &      'Error while executing a scalar function',fun,' ')
622    ENDIF
623  ENDIF
624!------------------------
625END SUBROUTINE mathop_r11
626!-
627!=== FUNCTIONS (only one argument)
628!-
629INTEGER FUNCTION ma_sin_r11(nb,x,nbo,y)
630!---------------------------------------------------------------------
631  IMPLICIT NONE
632!-
633  INTEGER :: nb,nbo,i
634  REAL :: x(nb),y(nbo)
635!---------------------------------------------------------------------
636  DO i=1,nb
637    y(i) = SIN(x(i))
638  ENDDO
639!-
640  nbo = nb
641  ma_sin_r11 = 0
642!----------------------
643END FUNCTION ma_sin_r11
644!===
645INTEGER FUNCTION ma_cos_r11(nb,x,nbo,y)
646!---------------------------------------------------------------------
647  IMPLICIT NONE
648!-
649  INTEGER :: nb,nbo,i
650  REAL :: x(nb),y(nbo)
651!---------------------------------------------------------------------
652  DO i=1,nb
653    y(i) = COS(x(i))
654  ENDDO
655!-
656  nbo = nb
657  ma_cos_r11 = 0
658!----------------------
659END FUNCTION ma_cos_r11
660!===
661INTEGER FUNCTION ma_tan_r11(nb,x,nbo,y)
662!---------------------------------------------------------------------
663  IMPLICIT NONE
664!-
665  INTEGER :: nb,nbo,i
666  REAL :: x(nb),y(nbo)
667!---------------------------------------------------------------------
668  DO i=1,nb
669    y(i) = TAN(x(i))
670  ENDDO
671!-
672  nbo = nb
673  ma_tan_r11 = 0
674!----------------------
675END FUNCTION ma_tan_r11
676!===
677INTEGER FUNCTION ma_asin_r11(nb,x,nbo,y)
678!---------------------------------------------------------------------
679  IMPLICIT NONE
680!-
681  INTEGER :: nb,nbo,i
682  REAL :: x(nb),y(nbo)
683!---------------------------------------------------------------------
684  DO i=1,nb
685    y(i) = ASIN(x(i))
686  ENDDO
687!-
688  nbo = nb
689  ma_asin_r11 = 0
690!-----------------------
691END FUNCTION ma_asin_r11
692!===
693INTEGER FUNCTION ma_acos_r11(nb,x,nbo,y)
694!---------------------------------------------------------------------
695  IMPLICIT NONE
696!-
697  INTEGER :: nb,nbo,i
698  REAL :: x(nb),y(nbo)
699!---------------------------------------------------------------------
700  DO i=1,nb
701    y(i) = ACOS(x(i))
702  ENDDO
703!-
704  nbo = nb
705  ma_acos_r11 = 0
706!-----------------------
707END FUNCTION ma_acos_r11
708!===
709INTEGER FUNCTION ma_atan_r11(nb,x,nbo,y)
710!---------------------------------------------------------------------
711  IMPLICIT NONE
712!-
713  INTEGER :: nb,nbo,i
714  REAL :: x(nb),y(nbo)
715!---------------------------------------------------------------------
716  DO i=1,nb
717    y(i) = ATAN(x(i))
718  ENDDO
719!-
720  nbo = nb
721  ma_atan_r11 = 0
722!-----------------------
723END FUNCTION ma_atan_r11
724!===
725INTEGER FUNCTION ma_exp_r11(nb,x,nbo,y)
726!---------------------------------------------------------------------
727  IMPLICIT NONE
728!-
729  INTEGER :: nb,nbo,i
730  REAL :: x(nb),y(nbo)
731!---------------------------------------------------------------------
732  DO i=1,nb
733    y(i) = EXP(x(i))
734  ENDDO
735!-
736  nbo = nb
737  ma_exp_r11 = 0
738!----------------------
739END FUNCTION ma_exp_r11
740!===
741INTEGER FUNCTION ma_log_r11(nb,x,nbo,y)
742!---------------------------------------------------------------------
743  IMPLICIT NONE
744!-
745  INTEGER :: nb,nbo,i
746  REAL :: x(nb),y(nbo)
747!---------------------------------------------------------------------
748  DO i=1,nb
749    y(i) = log(x(i))
750  ENDDO
751!-
752  nbo = nb
753  ma_log_r11 = 0
754!----------------------
755END FUNCTION ma_log_r11
756!===
757INTEGER FUNCTION ma_sqrt_r11(nb,x,nbo,y)
758!---------------------------------------------------------------------
759  IMPLICIT NONE
760!-
761  INTEGER :: nb,nbo,i
762  REAL :: x(nb),y(nbo)
763!---------------------------------------------------------------------
764  DO i=1,nb
765    y(i) = SQRT(x(i))
766  ENDDO
767!-
768  nbo = nb
769  ma_sqrt_r11 = 0
770!-----------------------
771END FUNCTION ma_sqrt_r11
772!===
773INTEGER FUNCTION ma_abs_r11(nb,x,nbo,y)
774!---------------------------------------------------------------------
775  IMPLICIT NONE
776!-
777  INTEGER :: nb,nbo,i
778  REAL :: x(nb),y(nbo)
779!---------------------------------------------------------------------
780  DO i=1,nb
781    y(i) = ABS(x(i))
782  ENDDO
783!-
784  nbo = nb
785  ma_abs_r11 = 0
786!----------------------
787END FUNCTION ma_abs_r11
788!===
789INTEGER FUNCTION ma_chs_r11(nb,x,nbo,y)
790!---------------------------------------------------------------------
791  IMPLICIT NONE
792!-
793  INTEGER :: nb,nbo,i
794  REAL :: x(nb),y(nbo)
795!---------------------------------------------------------------------
796  DO i=1,nb
797    y(i) = x(i)*(-1.)
798  ENDDO
799!-
800  nbo = nb
801  ma_chs_r11 = 0
802!----------------------
803END FUNCTION ma_chs_r11
804!===
805INTEGER FUNCTION ma_cels_r11(nb,x,nbo,y)
806!---------------------------------------------------------------------
807  IMPLICIT NONE
808!-
809  INTEGER :: nb,nbo,i
810  REAL :: x(nb),y(nbo)
811!---------------------------------------------------------------------
812  DO i=1,nb
813    y(i) = x(i)-273.15
814  ENDDO
815!-
816  nbo = nb
817  ma_cels_r11 = 0
818!-----------------------
819END FUNCTION ma_cels_r11
820!===
821INTEGER FUNCTION ma_kelv_r11(nb,x,nbo,y)
822!---------------------------------------------------------------------
823  IMPLICIT NONE
824!-
825  INTEGER :: nb,nbo,i
826  REAL :: x(nb),y(nbo)
827!---------------------------------------------------------------------
828  DO i=1,nb
829    y(i) = x(i)+273.15
830  ENDDO
831!-
832  nbo = nb
833  ma_kelv_r11 = 0
834!-----------------------
835END FUNCTION ma_kelv_r11
836!===
837INTEGER FUNCTION ma_deg_r11(nb,x,nbo,y)
838!---------------------------------------------------------------------
839  IMPLICIT NONE
840!-
841  INTEGER :: nb,nbo,i
842  REAL :: x(nb),y(nbo)
843!---------------------------------------------------------------------
844  DO i=1,nb
845    y(i) = x(i)*57.29577951
846  ENDDO
847!-
848  nbo = nb
849  ma_deg_r11 = 0
850!-----------------------
851END FUNCTION ma_deg_r11
852!===
853INTEGER FUNCTION ma_rad_r11(nb,x,nbo,y)
854!---------------------------------------------------------------------
855  IMPLICIT NONE
856!-
857  INTEGER :: nb,nbo,i
858  REAL :: x(nb),y(nbo)
859!---------------------------------------------------------------------
860  DO i=1,nb
861    y(i) = x(i)*0.01745329252
862  ENDDO
863!-
864  nbo = nb
865  ma_rad_r11 = 0
866!----------------------
867END FUNCTION ma_rad_r11
868!===
869INTEGER FUNCTION ma_ident_r11(nb,x,nbo,y)
870!---------------------------------------------------------------------
871  IMPLICIT NONE
872!-
873  INTEGER :: nb,nbo,i
874  REAL :: x(nb),y(nbo)
875!---------------------------------------------------------------------
876  DO i=1,nb
877    y(i) = x(i)
878  ENDDO
879!-
880  nbo = nb
881  ma_ident_r11 = 0
882!------------------------
883END FUNCTION ma_ident_r11
884!-
885!=== OPERATIONS (two argument)
886!-
887INTEGER FUNCTION ma_add_r11(nb,x,s,nbo,y)
888!---------------------------------------------------------------------
889  IMPLICIT NONE
890!-
891  INTEGER :: nb,nbo
892  REAL :: x(nb),s,y(nbo)
893!-
894  INTEGER :: i
895!---------------------------------------------------------------------
896  DO i=1,nb
897    y(i) = x(i)+s
898  ENDDO
899!-
900  nbo = nb
901  ma_add_r11 = 0
902!-----------------------
903  END FUNCTION ma_add_r11
904!===
905INTEGER FUNCTION ma_sub_r11(nb,x,s,nbo,y)
906!---------------------------------------------------------------------
907  IMPLICIT NONE
908!-
909  INTEGER :: nb,nbo
910  REAL :: x(nb),s,y(nbo)
911!-
912  INTEGER :: i
913!---------------------------------------------------------------------
914  DO i=1,nb
915    y(i) = x(i)-s
916  ENDDO
917!-
918  nbo = nb
919  ma_sub_r11 = 0
920!----------------------
921END FUNCTION ma_sub_r11
922!===
923INTEGER FUNCTION ma_subi_r11(nb,x,s,nbo,y)
924!---------------------------------------------------------------------
925  IMPLICIT NONE
926!-
927  INTEGER :: nb,nbo
928  REAL :: x(nb),s,y(nbo)
929!-
930  INTEGER :: i
931!---------------------------------------------------------------------
932  DO i=1,nb
933    y(i) =  s-x(i)
934  ENDDO
935!-
936  nbo = nb
937  ma_subi_r11 = 0
938!-----------------------
939END FUNCTION ma_subi_r11
940!===
941INTEGER FUNCTION ma_mult_r11(nb,x,s,nbo,y)
942!---------------------------------------------------------------------
943  IMPLICIT NONE
944!-
945  INTEGER :: nb,nbo
946  REAL :: x(nb),s,y(nbo)
947!-
948  INTEGER :: i
949!---------------------------------------------------------------------
950  DO i=1,nb
951    y(i) = x(i)*s
952  ENDDO
953!-
954  nbo = nb
955  ma_mult_r11 = 0
956!-----------------------
957END FUNCTION ma_mult_r11
958!===
959INTEGER FUNCTION ma_div_r11(nb,x,s,nbo,y)
960!---------------------------------------------------------------------
961  IMPLICIT NONE
962!-
963  INTEGER :: nb,nbo
964  REAL :: x(nb),s,y(nbo)
965!-
966  INTEGER :: i
967!---------------------------------------------------------------------
968  DO i=1,nb
969    y(i) = x(i)/s
970  ENDDO
971!-
972  nbo = nb
973  ma_div_r11 = 0
974!-----------------------
975  END FUNCTION ma_div_r11
976!===
977INTEGER FUNCTION ma_divi_r11(nb,x,s,nbo,y)
978!---------------------------------------------------------------------
979  IMPLICIT NONE
980!-
981  INTEGER :: nb,nbo
982  REAL :: x(nb),s,y(nbo)
983!-
984  INTEGER :: i
985!---------------------------------------------------------------------
986  DO i=1,nb
987    y(i) = s/x(i)
988  ENDDO
989!-
990  nbo = nb
991  ma_divi_r11 = 0
992!-----------------------
993END FUNCTION ma_divi_r11
994!===
995INTEGER FUNCTION ma_power_r11(nb,x,s,nbo,y)
996!---------------------------------------------------------------------
997  IMPLICIT NONE
998!-
999  INTEGER :: nb,nbo
1000  REAL :: x(nb),s,y(nbo)
1001!-
1002  INTEGER :: i
1003!---------------------------------------------------------------------
1004  DO i=1,nb
1005    y(i) = x(i)**s
1006  ENDDO
1007!-
1008  nbo = nb
1009  ma_power_r11 = 0
1010!-----------------------
1011END FUNCTION ma_power_r11
1012!===
1013INTEGER FUNCTION ma_fumin_r11(nb,x,s,nbo,y)
1014!---------------------------------------------------------------------
1015  IMPLICIT NONE
1016!-
1017  INTEGER :: nb,nbo
1018  REAL :: x(nb),s,y(nbo)
1019!-
1020  INTEGER :: i
1021!---------------------------------------------------------------------
1022  DO i=1,nb
1023    y(i) = MIN(x(i),s)
1024  ENDDO
1025!-
1026  nbo = nb
1027  ma_fumin_r11 = 0
1028!------------------------
1029END FUNCTION ma_fumin_r11
1030!===
1031INTEGER FUNCTION ma_fumax_r11(nb,x,s,nbo,y)
1032!---------------------------------------------------------------------
1033  IMPLICIT NONE
1034!-
1035  INTEGER :: nb,nbo
1036  REAL :: x(nb),s,y(nbo)
1037!-
1038  INTEGER :: i
1039!---------------------------------------------------------------------
1040  DO i=1,nb
1041    y(i) = MAX(x(i),s)
1042  ENDDO
1043!-
1044  nbo = nb
1045  ma_fumax_r11 = 0
1046!------------------------
1047END FUNCTION ma_fumax_r11
1048!===
1049INTEGER FUNCTION ma_fuscat_r11(nb,x,nbi,ind,miss_val,nbo,y)
1050!---------------------------------------------------------------------
1051  IMPLICIT NONE
1052!-
1053  INTEGER :: nb,nbo,nbi
1054  INTEGER :: ind(nbi)
1055  REAL :: x(nb),miss_val,y(nbo)
1056!-
1057  INTEGER :: i,ii,ipos
1058!---------------------------------------------------------------------
1059  ma_fuscat_r11 = 0
1060!-
1061  y(1:nbo) = miss_val
1062!-
1063  IF (nbi <= nb) THEN
1064    ipos = 0
1065    DO i=1,nbi
1066      IF (ind(i) <= nbo .AND. ind(i) > 0) THEN
1067        ipos = ipos+1
1068        y(ind(i)) = x(ipos)
1069      ELSE
1070        IF (ind(i) > nbo) ma_fuscat_r11  = ma_fuscat_r11+1
1071      ENDIF
1072    ENDDO
1073!-- Repeat the data if needed
1074    IF (MINVAL(ind) < 0) THEN
1075      DO i=1,nbi
1076        IF (ind(i) <= 0) THEN
1077          DO ii=1,ABS(ind(i))-1
1078            IF (ind(i+1)+ii <= nbo) THEN
1079              y(ind(i+1)+ii) = y(ind(i+1))
1080            ELSE
1081              ma_fuscat_r11  = ma_fuscat_r11+1
1082            ENDIF
1083          ENDDO
1084        ENDIF
1085      ENDDO
1086    ENDIF
1087  ELSE
1088    ma_fuscat_r11  = 1
1089  ENDIF
1090!-------------------------
1091END FUNCTION ma_fuscat_r11
1092!===
1093INTEGER FUNCTION ma_fugath_r11(nb,x,nbi,ind,miss_val,nbo,y)
1094!---------------------------------------------------------------------
1095  IMPLICIT NONE
1096!-
1097  INTEGER :: nb,nbo,nbi
1098  INTEGER :: ind(nbi)
1099  REAL :: x(nb),miss_val,y(nbo)
1100!-
1101  INTEGER :: i,ipos
1102!---------------------------------------------------------------------
1103  IF (nbi <= nbo) THEN
1104    ma_fugath_r11 = 0
1105    y(1:nbo) = miss_val
1106    ipos = 0
1107    DO i=1,nbi
1108      IF (ipos+1 <= nbo) THEN
1109        IF (ind(i) > 0) THEN
1110          ipos = ipos+1
1111          y(ipos) = x(ind(i))
1112        ENDIF
1113      ELSE
1114        IF (ipos+1 > nbo) ma_fugath_r11  = ma_fugath_r11+1
1115      ENDIF
1116    ENDDO
1117  ELSE
1118    ma_fugath_r11 = 1
1119  ENDIF
1120!-
1121  nbo = ipos
1122!-------------------------
1123END FUNCTION ma_fugath_r11
1124!===
1125INTEGER FUNCTION ma_fufill_r11(nb,x,nbi,ind,miss_val,nbo,y)
1126!---------------------------------------------------------------------
1127  IMPLICIT NONE
1128!-
1129  INTEGER :: nb,nbo,nbi
1130  INTEGER :: ind(nbi)
1131  REAL :: x(nb),miss_val,y(nbo)
1132!-
1133  INTEGER :: i,ii,ipos
1134!---------------------------------------------------------------------
1135  ma_fufill_r11 = 0
1136!-
1137  IF (nbi <= nb) THEN
1138    ipos = 0
1139    DO i=1,nbi
1140      IF (ind(i) <= nbo .AND. ind(i) > 0) THEN
1141        ipos = ipos+1
1142        y(ind(i)) = x(ipos)
1143      ELSE
1144        IF (ind(i) > nbo) ma_fufill_r11  = ma_fufill_r11+1
1145      ENDIF
1146    ENDDO
1147!-- Repeat the data if needed
1148    IF (MINVAL(ind) < 0) THEN
1149      DO i=1,nbi
1150        IF (ind(i) <= 0) THEN
1151          DO ii=1,ABS(ind(i))-1
1152            IF (ind(i+1)+ii <= nbo) THEN
1153              y(ind(i+1)+ii) = y(ind(i+1))
1154            ELSE
1155              ma_fufill_r11  = ma_fufill_r11+1
1156            ENDIF
1157          ENDDO
1158        ENDIF
1159      ENDDO
1160    ENDIF
1161  ELSE
1162    ma_fufill_r11  = 1
1163  ENDIF
1164!-------------------------
1165END FUNCTION ma_fufill_r11
1166!===
1167INTEGER FUNCTION ma_fucoll_r11(nb,x,nbi,ind,miss_val,nbo,y)
1168!---------------------------------------------------------------------
1169  IMPLICIT NONE
1170!-
1171  INTEGER :: nb,nbo,nbi
1172  INTEGER :: ind(nbi)
1173  REAL :: x(nb),miss_val,y(nbo)
1174!-
1175  INTEGER :: i,ipos
1176!---------------------------------------------------------------------
1177  IF (nbi <= nbo) THEN
1178    ma_fucoll_r11 = 0
1179    ipos = 0
1180    DO i=1,nbi
1181      IF (ipos+1 <= nbo) THEN
1182        IF (ind(i) > 0) THEN
1183          ipos = ipos+1
1184          y(ipos) = x(ind(i))
1185        ENDIF
1186      ELSE
1187        IF (ipos+1 > nbo) ma_fucoll_r11  = ma_fucoll_r11+1
1188      ENDIF
1189    ENDDO
1190  ELSE
1191    ma_fucoll_r11 = 1
1192  ENDIF
1193!-
1194  nbo = ipos
1195!-------------------------
1196END FUNCTION ma_fucoll_r11
1197!===
1198INTEGER FUNCTION ma_fuundef_r11(nb,x,nbi,ind,miss_val,nbo,y)
1199!---------------------------------------------------------------------
1200  IMPLICIT NONE
1201!-
1202  INTEGER :: nb,nbo,nbi
1203  INTEGER :: ind(nbi)
1204  REAL :: x(nb),miss_val,y(nbo)
1205!-
1206  INTEGER :: i
1207!---------------------------------------------------------------------
1208  IF (nbi <= nbo .AND. nbo == nb) THEN
1209    ma_fuundef_r11 = 0
1210    DO i=1,nbo
1211      y(i) = x(i)
1212    ENDDO
1213    DO i=1,nbi
1214      IF (ind(i) <= nbo .AND. ind(i) > 0) THEN
1215        y(ind(i)) =  miss_val
1216      ELSE
1217        IF (ind(i) > nbo) ma_fuundef_r11  = ma_fuundef_r11+1
1218      ENDIF
1219    ENDDO
1220  ELSE
1221    ma_fuundef_r11 = 1
1222  ENDIF
1223!--------------------------
1224END FUNCTION ma_fuundef_r11
1225!===
1226INTEGER FUNCTION ma_fuonly_r11(nb,x,nbi,ind,miss_val,nbo,y)
1227!---------------------------------------------------------------------
1228  IMPLICIT NONE
1229!-
1230  INTEGER :: nb,nbo,nbi
1231  INTEGER :: ind(nbi)
1232  REAL :: x(nb),miss_val,y(nbo)
1233!-
1234  INTEGER :: i
1235!---------------------------------------------------------------------
1236  IF (     (nbi <= nbo).AND.(nbo == nb) &
1237 &    .AND.ALL(ind(1:nbi) <= nbo) ) THEN
1238    ma_fuonly_r11 = 0
1239    y(1:nbo) = miss_val
1240    DO i=1,nbi
1241      IF (ind(i) > 0) THEN
1242        y(ind(i)) =  x(ind(i))
1243      ENDIF
1244    ENDDO
1245  ELSE
1246    ma_fuonly_r11 = 1
1247  ENDIF
1248!-------------------------
1249END FUNCTION ma_fuonly_r11
1250!===
1251!===
1252SUBROUTINE mathop_r21 &
1253 &  (fun,nb,work_in,miss_val,nb_index,nindex,scal,nb_max,work_out)
1254!---------------------------------------------------------------------
1255!- This subroutines gives an interface to the various operations
1256!- which are allowed. The interface is general enough to allow its use
1257!- for other cases.
1258!-
1259!- INPUT
1260!-
1261!- fun      : function to be applied to the vector of data
1262!- nb       : Length of input vector
1263!- work_in  : Input vector of data (REAL)
1264!- miss_val : The value of the missing data flag (it has to be a
1265!-            maximum value, in f90 : huge( a real ))
1266!- nb_index : Length of index vector
1267!- nindex   : Vector of indices
1268!- scal     : A scalar value for vector/scalar operations
1269!- nb_max   : maximum length of output vector
1270!-
1271!- OUTPUT
1272!-
1273!- nb_max   : Actual length of output variable
1274!- work_out : Output vector after the operation was applied
1275!---------------------------------------------------------------------
1276  IMPLICIT NONE
1277!-
1278  CHARACTER(LEN=7) :: fun
1279  INTEGER :: nb(2),nb_max,nb_index
1280  INTEGER :: nindex(nb_index)
1281  REAL :: work_in(nb(1),nb(2)),scal,miss_val
1282  REAL :: work_out(nb_max)
1283!-
1284  INTEGER :: ierr
1285!-
1286  INTRINSIC SIN,COS,TAN,ASIN,ACOS,ATAN,EXP,LOG,SQRT,ABS
1287!---------------------------------------------------------------------
1288  ierr = 0
1289!-
1290  IF (scal >= miss_val-1.) THEN
1291    IF (INDEX(indexfu,fun(1:LEN_TRIM(fun))) == 0) THEN
1292      SELECT CASE (fun)
1293      CASE('sin')
1294        ierr = ma_sin_r21(nb,work_in,nb_max,work_out)
1295      CASE('cos')
1296        ierr = ma_cos_r21(nb,work_in,nb_max,work_out)
1297      CASE('tan')
1298        ierr = ma_tan_r21(nb,work_in,nb_max,work_out)
1299      CASE('asin')
1300        ierr = ma_asin_r21(nb,work_in,nb_max,work_out)
1301      CASE('acos')
1302        ierr = ma_acos_r21(nb,work_in,nb_max,work_out)
1303      CASE('atan')
1304        ierr = ma_atan_r21(nb,work_in,nb_max,work_out)
1305      CASE('exp')
1306        ierr = ma_exp_r21(nb,work_in,nb_max,work_out)
1307      CASE('log')
1308        ierr = ma_log_r21(nb,work_in,nb_max,work_out)
1309      CASE('sqrt')
1310        ierr = ma_sqrt_r21(nb,work_in,nb_max,work_out)
1311      CASE('chs')
1312        ierr = ma_chs_r21(nb,work_in,nb_max,work_out)
1313      CASE('abs')
1314        ierr = ma_abs_r21(nb,work_in,nb_max,work_out)
1315      CASE('cels')
1316        ierr = ma_cels_r21(nb,work_in,nb_max,work_out)
1317      CASE('kelv')
1318        ierr = ma_kelv_r21(nb,work_in,nb_max,work_out)
1319      CASE('deg')
1320        ierr = ma_deg_r21(nb,work_in,nb_max,work_out)
1321      CASE('rad')
1322        ierr = ma_rad_r21(nb,work_in,nb_max,work_out)
1323      CASE('ident')
1324        ierr = ma_ident_r21(nb,work_in,nb_max,work_out)
1325      CASE DEFAULT
1326        CALL ipslerr(3,"mathop", &
1327 &        'scalar variable undefined and no indexing', &
1328 &        'but still unknown function',fun)
1329      END SELECT
1330      IF (ierr > 0) THEN
1331        CALL ipslerr(3,"mathop", &
1332 &        'Error while executing a simple function',fun,' ')
1333      ENDIF
1334    ELSE
1335      SELECT CASE (fun)
1336      CASE('gather')
1337        ierr = ma_fugath_r21(nb,work_in,nb_index,nindex, &
1338 &                            miss_val,nb_max,work_out)
1339      CASE('scatter')
1340        IF (nb_index > (nb(1)*nb(2)) ) THEN
1341          work_out(1:nb_max) = miss_val
1342          ierr=1
1343        ELSE
1344          ierr = ma_fuscat_r21(nb,work_in,nb_index,nindex, &
1345 &                             miss_val,nb_max,work_out)
1346        ENDIF
1347      CASE('coll')
1348        ierr = ma_fucoll_r21(nb,work_in,nb_index,nindex, &
1349 &                           miss_val,nb_max,work_out)
1350      CASE('fill')
1351        ierr = ma_fufill_r21(nb,work_in,nb_index,nindex, &
1352 &                           miss_val,nb_max,work_out)
1353      CASE('undef')
1354        ierr = ma_fuundef_r21(nb,work_in,nb_index,nindex, &
1355 &                           miss_val,nb_max,work_out)
1356      CASE('only')
1357        ierr = ma_fuonly_r21(nb,work_in,nb_index,nindex, &
1358 &                           miss_val,nb_max,work_out)
1359      CASE DEFAULT
1360        CALL ipslerr(3,"mathop", &
1361 &        'scalar variable undefined and indexing', &
1362 &        'was requested but with unknown function',fun)
1363      END SELECT
1364      IF (ierr > 0) THEN
1365        CALL ipslerr(3,"mathop_r21", &
1366 &        'Error while executing an indexing function',fun,' ')
1367      ENDIF
1368    ENDIF
1369  ELSE
1370    SELECT CASE (fun)
1371    CASE('fumin')
1372      ierr = ma_fumin_r21(nb,work_in,scal,nb_max,work_out)
1373    CASE('fumax')
1374      ierr = ma_fumax_r21(nb,work_in,scal,nb_max,work_out)
1375    CASE('add')
1376      ierr = ma_add_r21(nb,work_in,scal,nb_max,work_out)
1377    CASE('subi')
1378      ierr = ma_subi_r21(nb,work_in,scal,nb_max,work_out)
1379    CASE('sub')
1380      ierr = ma_sub_r21(nb,work_in,scal,nb_max,work_out)
1381    CASE('mult')
1382      ierr = ma_mult_r21(nb,work_in,scal,nb_max,work_out)
1383    CASE('div')
1384      ierr = ma_div_r21(nb,work_in,scal,nb_max,work_out)
1385    CASE('divi')
1386      ierr = ma_divi_r21(nb,work_in,scal,nb_max,work_out)
1387    CASE('power')
1388      ierr = ma_power_r21(nb,work_in,scal,nb_max,work_out)
1389    CASE DEFAULT
1390      CALL ipslerr(3,"mathop", &
1391 &      'Unknown operation with a scalar',fun,' ')
1392    END SELECT
1393    IF (ierr > 0) THEN
1394      CALL ipslerr(3,"mathop", &
1395 &      'Error while executing a scalar function',fun,' ')
1396    ENDIF
1397  ENDIF
1398!------------------------
1399END SUBROUTINE mathop_r21
1400!-
1401!=== FUNCTIONS (only one argument)
1402!-
1403INTEGER FUNCTION ma_sin_r21(nb,x,nbo,y)
1404!---------------------------------------------------------------------
1405  IMPLICIT NONE
1406!-
1407  INTEGER :: nb(2),nbo,i,j,ij
1408  REAL :: x(nb(1),nb(2)),y(nbo)
1409!---------------------------------------------------------------------
1410  ij = 0
1411  DO j=1,nb(2)
1412    DO i=1,nb(1)
1413      ij = ij+1
1414      y(ij) = SIN(x(i,j))
1415    ENDDO
1416  ENDDO
1417!-
1418  nbo = nb(1)*nb(2)
1419  ma_sin_r21 = 0
1420!----------------------
1421END FUNCTION ma_sin_r21
1422!===
1423INTEGER FUNCTION ma_cos_r21(nb,x,nbo,y)
1424!---------------------------------------------------------------------
1425  IMPLICIT NONE
1426!-
1427  INTEGER :: nb(2),nbo,i,j,ij
1428  REAL :: x(nb(1),nb(2)),y(nbo)
1429!---------------------------------------------------------------------
1430  ij = 0
1431  DO j=1,nb(2)
1432    DO i=1,nb(1)
1433      ij = ij+1
1434      y(ij) = COS(x(i,j))
1435    ENDDO
1436  ENDDO
1437!-
1438  nbo = nb(1)*nb(2)
1439  ma_cos_r21 = 0
1440!----------------------
1441END FUNCTION ma_cos_r21
1442!===
1443INTEGER FUNCTION ma_tan_r21(nb,x,nbo,y)
1444!---------------------------------------------------------------------
1445  IMPLICIT NONE
1446!-
1447  INTEGER :: nb(2),nbo,i,j,ij
1448  REAL :: x(nb(1),nb(2)),y(nbo)
1449!---------------------------------------------------------------------
1450  ij = 0
1451  DO j=1,nb(2)
1452    DO i=1,nb(1)
1453      ij = ij+1
1454      y(ij) = TAN(x(i,j))
1455    ENDDO
1456  ENDDO
1457!-
1458  nbo = nb(1)*nb(2)
1459  ma_tan_r21 = 0
1460!----------------------
1461END FUNCTION ma_tan_r21
1462!===
1463  INTEGER FUNCTION ma_asin_r21(nb,x,nbo,y)
1464!---------------------------------------------------------------------
1465  IMPLICIT NONE
1466!-
1467  INTEGER :: nb(2),nbo,i,j,ij
1468  REAL :: x(nb(1),nb(2)),y(nbo)
1469!---------------------------------------------------------------------
1470  ij = 0
1471  DO j=1,nb(2)
1472    DO i=1,nb(1)
1473      ij = ij+1
1474      y(ij) = ASIN(x(i,j))
1475    ENDDO
1476  ENDDO
1477!-
1478  nbo = nb(1)*nb(2)
1479  ma_asin_r21 = 0
1480!-----------------------
1481END FUNCTION ma_asin_r21
1482!===
1483INTEGER FUNCTION ma_acos_r21(nb,x,nbo,y)
1484!---------------------------------------------------------------------
1485  IMPLICIT NONE
1486!-
1487  INTEGER :: nb(2),nbo,i,j,ij
1488  REAL :: x(nb(1),nb(2)),y(nbo)
1489!---------------------------------------------------------------------
1490  ij = 0
1491  DO j=1,nb(2)
1492    DO i=1,nb(1)
1493      ij = ij+1
1494      y(ij) = ACOS(x(i,j))
1495    ENDDO
1496  ENDDO
1497!-
1498  nbo = nb(1)*nb(2)
1499  ma_acos_r21 = 0
1500!-----------------------
1501END FUNCTION ma_acos_r21
1502!===
1503INTEGER FUNCTION ma_atan_r21(nb,x,nbo,y)
1504!---------------------------------------------------------------------
1505  IMPLICIT NONE
1506!-
1507  INTEGER :: nb(2),nbo,i,j,ij
1508  REAL :: x(nb(1),nb(2)),y(nbo)
1509!---------------------------------------------------------------------
1510  ij = 0
1511  DO j=1,nb(2)
1512    DO i=1,nb(1)
1513      ij = ij+1
1514      y(ij) = ATAN(x(i,j))
1515    ENDDO
1516  ENDDO
1517!-
1518  nbo = nb(1)*nb(2)
1519  ma_atan_r21 = 0
1520!-----------------------
1521END FUNCTION ma_atan_r21
1522!===
1523INTEGER FUNCTION ma_exp_r21(nb,x,nbo,y)
1524!---------------------------------------------------------------------
1525  IMPLICIT NONE
1526!-
1527  INTEGER :: nb(2),nbo,i,j,ij
1528  REAL :: x(nb(1),nb(2)),y(nbo)
1529!---------------------------------------------------------------------
1530  ij = 0
1531  DO j=1,nb(2)
1532    DO i=1,nb(1)
1533      ij = ij+1
1534      y(ij) = EXP(x(i,j))
1535    ENDDO
1536  ENDDO
1537!-
1538  nbo = nb(1)*nb(2)
1539  ma_exp_r21 = 0
1540!----------------------
1541END FUNCTION ma_exp_r21
1542!===
1543INTEGER FUNCTION ma_log_r21(nb,x,nbo,y)
1544!---------------------------------------------------------------------
1545  IMPLICIT NONE
1546!-
1547  INTEGER :: nb(2),nbo,i,j,ij
1548  REAL :: x(nb(1),nb(2)),y(nbo)
1549!---------------------------------------------------------------------
1550  ij = 0
1551  DO j=1,nb(2)
1552    DO i=1,nb(1)
1553      ij = ij+1
1554      y(ij) = LOG(x(i,j))
1555    ENDDO
1556  ENDDO
1557!-
1558  nbo = nb(1)*nb(2)
1559  ma_log_r21 = 0
1560!----------------------
1561END FUNCTION ma_log_r21
1562!===
1563INTEGER FUNCTION ma_sqrt_r21(nb,x,nbo,y)
1564!---------------------------------------------------------------------
1565  IMPLICIT NONE
1566!-
1567  INTEGER :: nb(2),nbo,i,j,ij
1568  REAL :: x(nb(1),nb(2)),y(nbo)
1569!---------------------------------------------------------------------
1570  ij = 0
1571  DO j=1,nb(2)
1572    DO i=1,nb(1)
1573      ij = ij+1
1574      y(ij) = SQRT(x(i,j))
1575    ENDDO
1576  ENDDO
1577!-
1578  nbo = nb(1)*nb(2)
1579  ma_sqrt_r21 = 0
1580!-----------------------
1581END FUNCTION ma_sqrt_r21
1582!===
1583INTEGER FUNCTION ma_abs_r21(nb,x,nbo,y)
1584!---------------------------------------------------------------------
1585  IMPLICIT NONE
1586!-
1587  INTEGER :: nb(2),nbo,i,j,ij
1588  REAL :: x(nb(1),nb(2)),y(nbo)
1589!---------------------------------------------------------------------
1590  ij = 0
1591  DO j=1,nb(2)
1592    DO i=1,nb(1)
1593      ij = ij+1
1594      y(ij) = ABS(x(i,j))
1595    ENDDO
1596  ENDDO
1597!-
1598  nbo = nb(1)*nb(2)
1599  ma_abs_r21 = 0
1600!----------------------
1601END FUNCTION ma_abs_r21
1602!===
1603INTEGER FUNCTION ma_chs_r21(nb,x,nbo,y)
1604!---------------------------------------------------------------------
1605  IMPLICIT NONE
1606!-
1607  INTEGER :: nb(2),nbo,i,j,ij
1608  REAL :: x(nb(1),nb(2)),y(nbo)
1609!---------------------------------------------------------------------
1610  ij = 0
1611  DO j=1,nb(2)
1612    DO i=1,nb(1)
1613      ij = ij+1
1614      y(ij) = x(i,j)*(-1.)
1615    ENDDO
1616  ENDDO
1617!-
1618  nbo = nb(1)*nb(2)
1619  ma_chs_r21 = 0
1620!----------------------
1621END FUNCTION ma_chs_r21
1622!===
1623INTEGER FUNCTION ma_cels_r21(nb,x,nbo,y)
1624!---------------------------------------------------------------------
1625  IMPLICIT NONE
1626!-
1627  INTEGER :: nb(2),nbo,i,j,ij
1628  REAL :: x(nb(1),nb(2)),y(nbo)
1629!---------------------------------------------------------------------
1630  ij = 0
1631  DO j=1,nb(2)
1632    DO i=1,nb(1)
1633      ij = ij+1
1634      y(ij) = x(i,j)-273.15
1635    ENDDO
1636  ENDDO
1637!-
1638  nbo = nb(1)*nb(2)
1639  ma_cels_r21 = 0
1640!-----------------------
1641END FUNCTION ma_cels_r21
1642!===
1643INTEGER FUNCTION ma_kelv_r21(nb,x,nbo,y)
1644!---------------------------------------------------------------------
1645  IMPLICIT NONE
1646!-
1647  INTEGER :: nb(2),nbo,i,j,ij
1648  REAL :: x(nb(1),nb(2)),y(nbo)
1649!---------------------------------------------------------------------
1650  ij = 0
1651  DO j=1,nb(2)
1652    DO i=1,nb(1)
1653      ij = ij+1
1654      y(ij) = x(i,j)+273.15
1655    ENDDO
1656  ENDDO
1657!-
1658  nbo = nb(1)*nb(2)
1659  ma_kelv_r21 = 0
1660!-----------------------
1661END FUNCTION ma_kelv_r21
1662!===
1663INTEGER FUNCTION ma_deg_r21(nb,x,nbo,y)
1664!---------------------------------------------------------------------
1665  IMPLICIT NONE
1666!-
1667  INTEGER :: nb(2),nbo,i,j,ij
1668  REAL :: x(nb(1),nb(2)),y(nbo)
1669!---------------------------------------------------------------------
1670  ij = 0
1671  DO j=1,nb(2)
1672    DO i=1,nb(1)
1673      ij = ij+1
1674      y(ij) = x(i,j)*57.29577951
1675    ENDDO
1676  ENDDO
1677!-
1678  nbo = nb(1)*nb(2)
1679  ma_deg_r21 = 0
1680!----------------------
1681END FUNCTION ma_deg_r21
1682!===
1683INTEGER FUNCTION ma_rad_r21(nb,x,nbo,y)
1684!---------------------------------------------------------------------
1685  IMPLICIT NONE
1686!-
1687  INTEGER :: nb(2),nbo,i,j,ij
1688  REAL :: x(nb(1),nb(2)),y(nbo)
1689!---------------------------------------------------------------------
1690  ij = 0
1691  DO j=1,nb(2)
1692    DO i=1,nb(1)
1693      ij = ij+1
1694      y(ij) = x(i,j)*0.01745329252
1695    ENDDO
1696  ENDDO
1697!-
1698  nbo = nb(1)*nb(2)
1699  ma_rad_r21 = 0
1700!----------------------
1701END FUNCTION ma_rad_r21
1702!===
1703INTEGER FUNCTION ma_ident_r21(nb,x,nbo,y)
1704!---------------------------------------------------------------------
1705  IMPLICIT NONE
1706!-
1707  INTEGER :: nb(2),nbo,i,j,ij
1708  REAL :: x(nb(1),nb(2)),y(nbo)
1709!---------------------------------------------------------------------
1710  ij = 0
1711  DO j=1,nb(2)
1712    DO i=1,nb(1)
1713      ij = ij+1
1714      y(ij) = x(i,j)
1715    ENDDO
1716  ENDDO
1717!-
1718  nbo = nb(1)*nb(2)
1719  ma_ident_r21 = 0
1720!------------------------
1721END FUNCTION ma_ident_r21
1722!-
1723!=== OPERATIONS (two argument)
1724!-
1725INTEGER FUNCTION ma_add_r21(nb,x,s,nbo,y)
1726!---------------------------------------------------------------------
1727  IMPLICIT NONE
1728!-
1729  INTEGER :: nb(2),nbo
1730  REAL :: x(nb(1),nb(2)),s,y(nbo)
1731!-
1732  INTEGER :: i,j,ij
1733!---------------------------------------------------------------------
1734  ij = 0
1735  DO j=1,nb(2)
1736    DO i=1,nb(1)
1737      ij = ij+1
1738      y(ij) = x(i,j)+s
1739    ENDDO
1740  ENDDO
1741!-
1742  nbo = nb(1)*nb(2)
1743  ma_add_r21 = 0
1744!----------------------
1745END FUNCTION ma_add_r21
1746!===
1747INTEGER FUNCTION ma_sub_r21(nb,x,s,nbo,y)
1748!---------------------------------------------------------------------
1749  IMPLICIT NONE
1750!-
1751  INTEGER :: nb(2),nbo
1752  REAL :: x(nb(1),nb(2)),s,y(nbo)
1753!-
1754  INTEGER :: i,j,ij
1755!---------------------------------------------------------------------
1756  ij = 0
1757  DO j=1,nb(2)
1758    DO i=1,nb(1)
1759      ij = ij+1
1760      y(ij) = x(i,j)-s
1761    ENDDO
1762  ENDDO
1763!-
1764  nbo = nb(1)*nb(2)
1765  ma_sub_r21 = 0
1766!----------------------
1767END FUNCTION ma_sub_r21
1768!===
1769INTEGER FUNCTION ma_subi_r21(nb,x,s,nbo,y)
1770!---------------------------------------------------------------------
1771  IMPLICIT NONE
1772!-
1773  INTEGER :: nb(2),nbo
1774  REAL :: x(nb(1),nb(2)),s,y(nbo)
1775!-
1776  INTEGER :: i,j,ij
1777!---------------------------------------------------------------------
1778  ij = 0
1779  DO j=1,nb(2)
1780    DO i=1,nb(1)
1781      ij = ij+1
1782      y(ij) =  s-x(i,j)
1783    ENDDO
1784  ENDDO
1785!-
1786  nbo = nb(1)*nb(2)
1787  ma_subi_r21 = 0
1788!-----------------------
1789END FUNCTION ma_subi_r21
1790!===
1791INTEGER FUNCTION ma_mult_r21(nb,x,s,nbo,y)
1792!---------------------------------------------------------------------
1793  IMPLICIT NONE
1794!-
1795  INTEGER :: nb(2),nbo
1796  REAL :: x(nb(1),nb(2)),s,y(nbo)
1797!-
1798  INTEGER :: i,j,ij
1799!---------------------------------------------------------------------
1800  ij = 0
1801  DO j=1,nb(2)
1802    DO i=1,nb(1)
1803      ij = ij+1
1804      y(ij) = x(i,j)*s
1805    ENDDO
1806  ENDDO
1807!-
1808  nbo = nb(1)*nb(2)
1809  ma_mult_r21 = 0
1810!-----------------------
1811END FUNCTION ma_mult_r21
1812!===
1813INTEGER FUNCTION ma_div_r21(nb,x,s,nbo,y)
1814!---------------------------------------------------------------------
1815  IMPLICIT NONE
1816!-
1817  INTEGER :: nb(2),nbo
1818  REAL :: x(nb(1),nb(2)),s,y(nbo)
1819!-
1820  INTEGER :: i,j,ij
1821!---------------------------------------------------------------------
1822  ij = 0
1823  DO j=1,nb(2)
1824    DO i=1,nb(1)
1825      ij = ij+1
1826      y(ij) = x(i,j)/s
1827    ENDDO
1828  ENDDO
1829!-
1830  nbo = nb(1)*nb(2)
1831  ma_div_r21 = 0
1832!----------------------
1833END FUNCTION ma_div_r21
1834!===
1835INTEGER FUNCTION ma_divi_r21(nb,x,s,nbo,y)
1836!---------------------------------------------------------------------
1837  IMPLICIT NONE
1838!-
1839  INTEGER :: nb(2),nbo
1840  REAL :: x(nb(1),nb(2)),s,y(nbo)
1841!-
1842  INTEGER :: i,j,ij
1843!---------------------------------------------------------------------
1844  ij = 0
1845  DO j=1,nb(2)
1846    DO i=1,nb(1)
1847      ij = ij+1
1848      y(ij) = s/x(i,j)
1849    ENDDO
1850  ENDDO
1851!-
1852  nbo = nb(1)*nb(2)
1853  ma_divi_r21 = 0
1854!-----------------------
1855END FUNCTION ma_divi_r21
1856!===
1857INTEGER FUNCTION ma_power_r21(nb,x,s,nbo,y)
1858!---------------------------------------------------------------------
1859  IMPLICIT NONE
1860!-
1861  INTEGER :: nb(2),nbo
1862  REAL :: x(nb(1),nb(2)),s,y(nbo)
1863!-
1864  INTEGER :: i,j,ij
1865!---------------------------------------------------------------------
1866  ij = 0
1867  DO j=1,nb(2)
1868    DO i=1,nb(1)
1869      ij = ij+1
1870      y(ij) = x(i,j) ** s
1871    ENDDO
1872  ENDDO
1873!-
1874  nbo = nb(1)*nb(2)
1875  ma_power_r21 = 0
1876!------------------------
1877END FUNCTION ma_power_r21
1878!===
1879INTEGER FUNCTION ma_fumin_r21(nb,x,s,nbo,y)
1880!---------------------------------------------------------------------
1881  IMPLICIT NONE
1882!-
1883  INTEGER :: nb(2),nbo
1884  REAL :: x(nb(1),nb(2)),s,y(nbo)
1885!-
1886  INTEGER :: i,j,ij
1887!---------------------------------------------------------------------
1888  ij = 0
1889  DO j=1,nb(2)
1890    DO i=1,nb(1)
1891      ij = ij+1
1892      y(ij) = MIN(x(i,j),s)
1893    ENDDO
1894  ENDDO
1895!-
1896  nbo = nb(1)*nb(2)
1897  ma_fumin_r21 = 0
1898!------------------------
1899END FUNCTION ma_fumin_r21
1900!===
1901INTEGER FUNCTION ma_fumax_r21(nb,x,s,nbo,y)
1902!---------------------------------------------------------------------
1903  IMPLICIT NONE
1904!-
1905  INTEGER :: nb(2),nbo
1906  REAL :: x(nb(1),nb(2)),s,y(nbo)
1907!-
1908  INTEGER :: i,j,ij
1909!---------------------------------------------------------------------
1910  ij = 0
1911  DO j=1,nb(2)
1912    DO i=1,nb(1)
1913      ij = ij+1
1914      y(ij) = MAX(x(i,j),s)
1915    ENDDO
1916  ENDDO
1917!-
1918  nbo = nb(1)*nb(2)
1919  ma_fumax_r21 = 0
1920!------------------------
1921END FUNCTION ma_fumax_r21
1922!===
1923INTEGER FUNCTION ma_fuscat_r21(nb,x,nbi,ind,miss_val,nbo,y)
1924!---------------------------------------------------------------------
1925  IMPLICIT NONE
1926!-
1927  INTEGER :: nb(2),nbo,nbi
1928  INTEGER :: ind(nbi)
1929  REAL :: x(nb(1),nb(2)),miss_val,y(nbo)
1930!-
1931  INTEGER :: i,j,ij,ii,ipos
1932!---------------------------------------------------------------------
1933  ma_fuscat_r21 = 0
1934!-
1935  y(1:nbo) = miss_val
1936!-
1937  IF (nbi <= nb(1)*nb(2)) THEN
1938    ipos = 0
1939    DO ij=1,nbi
1940       IF (ind(ij) <= nbo .AND. ind(ij) > 0) THEN
1941         ipos = ipos+1
1942         j = ((ipos-1)/nb(1))+1
1943         i = (ipos-(j-1)*nb(1))
1944         y(ind(ij)) = x(i,j)
1945       ELSE
1946         IF (ind(ij) > nbo) ma_fuscat_r21  = ma_fuscat_r21+1
1947       ENDIF
1948    ENDDO
1949!-- Repeat the data if needed
1950    IF (MINVAL(ind) < 0) THEN
1951      DO i=1,nbi
1952        IF (ind(i) <= 0) THEN
1953          DO ii=1,ABS(ind(i))-1
1954            IF (ind(i+1)+ii <= nbo) THEN
1955              y(ind(i+1)+ii) = y(ind(i+1))
1956            ELSE
1957              ma_fuscat_r21  = ma_fuscat_r21+1
1958            ENDIF
1959          ENDDO
1960        ENDIF
1961      ENDDO
1962    ENDIF
1963  ELSE
1964    ma_fuscat_r21  = 1
1965  ENDIF
1966!-------------------------
1967END FUNCTION ma_fuscat_r21
1968!===
1969INTEGER FUNCTION ma_fugath_r21(nb,x,nbi,ind,miss_val,nbo,y)
1970!---------------------------------------------------------------------
1971  IMPLICIT NONE
1972!-
1973  INTEGER :: nb(2),nbo,nbi
1974  INTEGER :: ind(nbi)
1975  REAL :: x(nb(1),nb(2)),miss_val,y(nbo)
1976!-
1977  INTEGER :: i,j,ij,ipos
1978!---------------------------------------------------------------------
1979  IF (nbi <= nbo) THEN
1980    ma_fugath_r21 = 0
1981    y(1:nbo) = miss_val
1982    ipos = 0
1983    DO ij=1,nbi
1984      IF (ipos+1 <= nbo) THEN
1985        IF (ind(ij) > 0) THEN
1986          j = ((ind(ij)-1)/nb(1))+1
1987          i = (ind(ij)-(j-1)*nb(1))
1988          ipos = ipos+1
1989          y(ipos) = x(i,j)
1990        ENDIF
1991      ELSE
1992        IF (ipos+1 > nbo) ma_fugath_r21  = ma_fugath_r21+1
1993      ENDIF
1994    ENDDO
1995  ELSE
1996    ma_fugath_r21 = 1
1997  ENDIF
1998  nbo = ipos
1999!-------------------------
2000END FUNCTION ma_fugath_r21
2001!===
2002INTEGER FUNCTION ma_fufill_r21(nb,x,nbi,ind,miss_val,nbo,y)
2003!---------------------------------------------------------------------
2004  IMPLICIT NONE
2005!-
2006  INTEGER :: nb(2),nbo,nbi
2007  INTEGER :: ind(nbi)
2008  REAL :: x(nb(1),nb(2)),miss_val,y(nbo)
2009!-
2010  INTEGER :: i,j,ij,ii,ipos
2011!---------------------------------------------------------------------
2012  ma_fufill_r21 = 0
2013!-
2014  IF (nbi <= nb(1)*nb(2)) THEN
2015    ipos = 0
2016    DO ij=1,nbi
2017      IF (ind(ij) <= nbo .AND. ind(ij) > 0) THEN
2018        ipos = ipos+1
2019        j = ((ipos-1)/nb(1))+1
2020        i = (ipos-(j-1)*nb(1))
2021        y(ind(ij)) = x(i,j)
2022      ELSE
2023        IF (ind(ij) > nbo) ma_fufill_r21  = ma_fufill_r21+1
2024      ENDIF
2025    ENDDO
2026!-- Repeat the data if needed
2027    IF (MINVAL(ind) < 0) THEN
2028      DO i=1,nbi
2029        IF (ind(i) <= 0) THEN
2030          DO ii=1,ABS(ind(i))-1
2031            IF (ind(i+1)+ii <= nbo) THEN
2032              y(ind(i+1)+ii) = y(ind(i+1))
2033            ELSE
2034              ma_fufill_r21  = ma_fufill_r21+1
2035            ENDIF
2036          ENDDO
2037        ENDIF
2038      ENDDO
2039    ENDIF
2040  ELSE
2041    ma_fufill_r21  = 1
2042  ENDIF
2043!-------------------------
2044END FUNCTION ma_fufill_r21
2045!===
2046INTEGER FUNCTION ma_fucoll_r21(nb,x,nbi,ind,miss_val,nbo,y)
2047!---------------------------------------------------------------------
2048  IMPLICIT NONE
2049!-
2050  INTEGER :: nb(2),nbo,nbi
2051  INTEGER :: ind(nbi)
2052  REAL :: x(nb(1),nb(2)),miss_val,y(nbo)
2053!-
2054  INTEGER :: i,j,ij,ipos
2055!---------------------------------------------------------------------
2056  IF (nbi <= nbo) THEN
2057    ma_fucoll_r21 = 0
2058    ipos = 0
2059    DO ij=1,nbi
2060      IF (ipos+1 <= nbo) THEN
2061        IF (ind(ij) > 0) THEN
2062          j = ((ind(ij)-1)/nb(1))+1
2063          i = (ind(ij)-(j-1)*nb(1))
2064          ipos = ipos+1
2065          y(ipos) = x(i,j)
2066        ENDIF
2067      ELSE
2068        IF (ipos+1 > nbo) ma_fucoll_r21  = ma_fucoll_r21+1
2069      ENDIF
2070    ENDDO
2071  ELSE
2072    ma_fucoll_r21 = 1
2073  ENDIF
2074  nbo = ipos
2075!-------------------------
2076END FUNCTION ma_fucoll_r21
2077!===
2078INTEGER FUNCTION ma_fuundef_r21(nb,x,nbi,ind,miss_val,nbo,y)
2079!---------------------------------------------------------------------
2080  IMPLICIT NONE
2081!-
2082  INTEGER :: nb(2),nbo,nbi
2083  INTEGER :: ind(nbi)
2084  REAL :: x(nb(1),nb(2)),miss_val,y(nbo)
2085!-
2086  INTEGER :: i,j,ij
2087!---------------------------------------------------------------------
2088  IF (nbi <= nbo .AND. nbo == nb(1)*nb(2)) THEN
2089    ma_fuundef_r21 = 0
2090    DO ij=1,nbo
2091      j = ((ij-1)/nb(1))+1
2092      i = (ij-(j-1)*nb(1))
2093      y(ij) = x(i,j)
2094    ENDDO
2095    DO i=1,nbi
2096      IF (ind(i) <= nbo .AND. ind(i) > 0) THEN
2097        y(ind(i)) =  miss_val
2098      ELSE
2099        IF (ind(i) > nbo) ma_fuundef_r21 = ma_fuundef_r21+1
2100      ENDIF
2101    ENDDO
2102  ELSE
2103    ma_fuundef_r21 = 1
2104  ENDIF
2105!--------------------------
2106END FUNCTION ma_fuundef_r21
2107!===
2108INTEGER FUNCTION ma_fuonly_r21(nb,x,nbi,ind,miss_val,nbo,y)
2109!---------------------------------------------------------------------
2110  IMPLICIT NONE
2111!-
2112  INTEGER :: nb(2),nbo,nbi
2113  INTEGER :: ind(nbi)
2114  REAL :: x(nb(1),nb(2)),miss_val,y(nbo)
2115!-
2116  INTEGER :: i,j,ij
2117!---------------------------------------------------------------------
2118  IF (     (nbi <= nbo).AND.(nbo == nb(1)*nb(2)) &
2119 &    .AND.ALL(ind(1:nbi) <= nbo) ) THEN
2120    ma_fuonly_r21 = 0
2121    y(1:nbo) = miss_val
2122    DO ij=1,nbi
2123      IF (ind(ij) > 0) THEN
2124        j = ((ind(ij)-1)/nb(1))+1
2125        i = (ind(ij)-(j-1)*nb(1))
2126        y(ind(ij)) =  x(i,j)
2127      ENDIF
2128    ENDDO
2129  ELSE
2130    ma_fuonly_r21 = 1
2131  ENDIF
2132!-------------------------
2133END FUNCTION ma_fuonly_r21
2134!===
2135!===
2136SUBROUTINE mathop_r31 &
2137 &  (fun,nb,work_in,miss_val,nb_index,nindex,scal,nb_max,work_out)
2138!---------------------------------------------------------------------
2139!- This subroutines gives an interface to the various operations
2140!- which are allowed. The interface is general enough to allow its use
2141!- for other cases.
2142!-
2143!- INPUT
2144!-
2145!- fun      : function to be applied to the vector of data
2146!- nb       : Length of input vector
2147!- work_in  : Input vector of data (REAL)
2148!- miss_val : The value of the missing data flag (it has to be a
2149!-            maximum value, in f90 : huge( a real ))
2150!- nb_index : Length of index vector
2151!- nindex   : Vector of indices
2152!- scal     : A scalar value for vector/scalar operations
2153!- nb_max   : maximum length of output vector
2154!-
2155!- OUTPUT
2156!-
2157!- nb_max   : Actual length of output variable
2158!- work_out : Output vector after the operation was applied
2159!---------------------------------------------------------------------
2160  IMPLICIT NONE
2161!-
2162  CHARACTER(LEN=7) :: fun
2163  INTEGER :: nb(3),nb_max,nb_index
2164  INTEGER :: nindex(nb_index)
2165  REAL :: work_in(nb(1),nb(2),nb(3)),scal,miss_val
2166  REAL :: work_out(nb_max)
2167!-
2168  INTEGER :: ierr
2169!-
2170  INTRINSIC SIN,COS,TAN,ASIN,ACOS,ATAN,EXP,LOG,SQRT,ABS
2171!---------------------------------------------------------------------
2172  ierr = 0
2173!-
2174  IF (scal >= miss_val-1.) THEN
2175    IF (INDEX(indexfu,fun(1:LEN_TRIM(fun))) == 0) THEN
2176      SELECT CASE (fun)
2177      CASE('sin')
2178        ierr = ma_sin_r31(nb,work_in,nb_max,work_out)
2179      CASE('cos')
2180        ierr = ma_cos_r31(nb,work_in,nb_max,work_out)
2181      CASE('tan')
2182        ierr = ma_tan_r31(nb,work_in,nb_max,work_out)
2183      CASE('asin')
2184        ierr = ma_asin_r31(nb,work_in,nb_max,work_out)
2185      CASE('acos')
2186        ierr = ma_acos_r31(nb,work_in,nb_max,work_out)
2187      CASE('atan')
2188        ierr = ma_atan_r31(nb,work_in,nb_max,work_out)
2189      CASE('exp')
2190        ierr = ma_exp_r31(nb,work_in,nb_max,work_out)
2191      CASE('log')
2192        ierr = ma_log_r31(nb,work_in,nb_max,work_out)
2193      CASE('sqrt')
2194        ierr = ma_sqrt_r31(nb,work_in,nb_max,work_out)
2195      CASE('chs')
2196        ierr = ma_chs_r31(nb,work_in,nb_max,work_out)
2197      CASE('abs')
2198        ierr = ma_abs_r31(nb,work_in,nb_max,work_out)
2199      CASE('cels')
2200        ierr = ma_cels_r31(nb,work_in,nb_max,work_out)
2201      CASE('kelv')
2202        ierr = ma_kelv_r31(nb,work_in,nb_max,work_out)
2203      CASE('deg')
2204        ierr = ma_deg_r31(nb,work_in,nb_max,work_out)
2205      CASE('rad')
2206        ierr = ma_rad_r31(nb,work_in,nb_max,work_out)
2207      CASE('ident')
2208        ierr = ma_ident_r31(nb,work_in,nb_max,work_out)
2209      CASE DEFAULT
2210        CALL ipslerr(3,"mathop", &
2211 &        'scalar variable undefined and no indexing', &
2212 &        'but still unknown function',fun)
2213      END SELECT
2214      IF (ierr > 0) THEN
2215         CALL ipslerr(3,"mathop", &
2216 &        'Error while executing a simple function',fun,' ')
2217      ENDIF
2218    ELSE
2219      SELECT CASE (fun)
2220      CASE('gather')
2221        ierr = ma_fugath_r31(nb,work_in,nb_index,nindex, &
2222 &                           miss_val,nb_max,work_out)
2223      CASE('scatter')
2224        IF (nb_index > (nb(1)*nb(2)*nb(3))) THEN
2225          work_out(1:nb_max) = miss_val
2226          ierr=1
2227        ELSE
2228          ierr = ma_fuscat_r31(nb,work_in,nb_index,nindex, &
2229 &                             miss_val,nb_max,work_out)
2230        ENDIF
2231      CASE('coll')
2232        ierr = ma_fucoll_r31(nb,work_in,nb_index,nindex, &
2233 &                           miss_val,nb_max,work_out)
2234      CASE('fill')
2235        ierr = ma_fufill_r31(nb,work_in,nb_index,nindex, &
2236 &                           miss_val,nb_max,work_out)
2237      CASE('undef')
2238        ierr = ma_fuundef_r31(nb,work_in,nb_index,nindex, &
2239 &                           miss_val,nb_max,work_out)
2240      CASE('only')
2241        ierr = ma_fuonly_r31(nb,work_in,nb_index,nindex, &
2242 &                           miss_val,nb_max,work_out)
2243      CASE DEFAULT
2244        CALL ipslerr(3,"mathop", &
2245 &        'scalar variable undefined and indexing', &
2246 &        'was requested but with unknown function',fun)
2247      END SELECT
2248      IF (ierr > 0) THEN
2249        CALL ipslerr(3,"mathop_r31", &
2250 &        'Error while executing an indexing function',fun,' ')
2251      ENDIF
2252    ENDIF
2253  ELSE
2254    SELECT CASE (fun)
2255    CASE('fumin')
2256      ierr = ma_fumin_r31(nb,work_in,scal,nb_max,work_out)
2257    CASE('fumax')
2258      ierr = ma_fumax_r31(nb,work_in,scal,nb_max,work_out)
2259    CASE('add')
2260      ierr = ma_add_r31(nb,work_in,scal,nb_max,work_out)
2261    CASE('subi')
2262      ierr = ma_subi_r31(nb,work_in,scal,nb_max,work_out)
2263    CASE('sub')
2264      ierr = ma_sub_r31(nb,work_in,scal,nb_max,work_out)
2265    CASE('mult')
2266      ierr = ma_mult_r31(nb,work_in,scal,nb_max,work_out)
2267    CASE('div')
2268      ierr = ma_div_r31(nb,work_in,scal,nb_max,work_out)
2269    CASE('divi')
2270      ierr = ma_divi_r31(nb,work_in,scal,nb_max,work_out)
2271    CASE('power')
2272      ierr = ma_power_r31(nb,work_in,scal,nb_max,work_out)
2273    CASE DEFAULT
2274      CALL ipslerr(3,"mathop", &
2275 &      'Unknown operation with a scalar',fun,' ')
2276    END SELECT
2277    IF (ierr > 0) THEN
2278      CALL ipslerr(3,"mathop", &
2279 &      'Error while executing a scalar function',fun,' ')
2280    ENDIF
2281  ENDIF
2282!------------------------
2283END SUBROUTINE mathop_r31
2284!-
2285!=== FUNCTIONS (only one argument)
2286!-
2287INTEGER FUNCTION ma_sin_r31(nb,x,nbo,y)
2288!---------------------------------------------------------------------
2289  IMPLICIT NONE
2290!-
2291  INTEGER :: nb(3),nbo,i,j,k,ij
2292  REAL :: x(nb(1),nb(2),nb(3)),y(nbo)
2293!---------------------------------------------------------------------
2294  ij = 0
2295  DO k=1,nb(3)
2296    DO j=1,nb(2)
2297      DO i=1,nb(1)
2298        ij = ij+1
2299        y(ij) = SIN(x(i,j,k))
2300      ENDDO
2301    ENDDO
2302  ENDDO
2303!-
2304  nbo = nb(1)*nb(2)*nb(3)
2305  ma_sin_r31 = 0
2306!----------------------
2307END FUNCTION ma_sin_r31
2308!===
2309INTEGER FUNCTION ma_cos_r31(nb,x,nbo,y)
2310!---------------------------------------------------------------------
2311  IMPLICIT NONE
2312!-
2313  INTEGER :: nb(3),nbo,i,j,k,ij
2314  REAL :: x(nb(1),nb(2),nb(3)),y(nbo)
2315!---------------------------------------------------------------------
2316  ij = 0
2317  DO k=1,nb(3)
2318    DO j=1,nb(2)
2319      DO i=1,nb(1)
2320        ij = ij+1
2321        y(ij) = COS(x(i,j,k))
2322      ENDDO
2323    ENDDO
2324  ENDDO
2325!-
2326  nbo = nb(1)*nb(2)*nb(3)
2327  ma_cos_r31 = 0
2328!----------------------
2329END FUNCTION ma_cos_r31
2330!===
2331INTEGER FUNCTION ma_tan_r31(nb,x,nbo,y)
2332!---------------------------------------------------------------------
2333  IMPLICIT NONE
2334!-
2335  INTEGER :: nb(3),nbo,i,j,k,ij
2336  REAL :: x(nb(1),nb(2),nb(3)),y(nbo)
2337!---------------------------------------------------------------------
2338  ij = 0
2339  DO k=1,nb(3)
2340    DO j=1,nb(2)
2341      DO i=1,nb(1)
2342        ij = ij+1
2343        y(ij) = TAN(x(i,j,k))
2344      ENDDO
2345    ENDDO
2346  ENDDO
2347!-
2348  nbo = nb(1)*nb(2)*nb(3)
2349  ma_tan_r31 = 0
2350!----------------------
2351END FUNCTION ma_tan_r31
2352!===
2353INTEGER FUNCTION ma_asin_r31(nb,x,nbo,y)
2354!---------------------------------------------------------------------
2355  IMPLICIT NONE
2356!-
2357  INTEGER :: nb(3),nbo,i,j,k,ij
2358  REAL :: x(nb(1),nb(2),nb(3)),y(nbo)
2359!---------------------------------------------------------------------
2360  ij = 0
2361  DO k=1,nb(3)
2362    DO j=1,nb(2)
2363      DO i=1,nb(1)
2364        ij = ij+1
2365        y(ij) = ASIN(x(i,j,k))
2366      ENDDO
2367    ENDDO
2368  ENDDO
2369!-
2370  nbo = nb(1)*nb(2)*nb(3)
2371  ma_asin_r31 = 0
2372!-----------------------
2373END FUNCTION ma_asin_r31
2374!===
2375INTEGER FUNCTION ma_acos_r31(nb,x,nbo,y)
2376!---------------------------------------------------------------------
2377  IMPLICIT NONE
2378!-
2379  INTEGER :: nb(3),nbo,i,j,k,ij
2380  REAL :: x(nb(1),nb(2),nb(3)),y(nbo)
2381!---------------------------------------------------------------------
2382  ij = 0
2383  DO k=1,nb(3)
2384    DO j=1,nb(2)
2385      DO i=1,nb(1)
2386        ij = ij+1
2387        y(ij) = ACOS(x(i,j,k))
2388      ENDDO
2389    ENDDO
2390  ENDDO
2391!-
2392  nbo = nb(1)*nb(2)*nb(3)
2393  ma_acos_r31 = 0
2394!-----------------------
2395END FUNCTION ma_acos_r31
2396!===
2397INTEGER FUNCTION ma_atan_r31(nb,x,nbo,y)
2398!---------------------------------------------------------------------
2399  IMPLICIT NONE
2400!-
2401  INTEGER :: nb(3),nbo,i,j,k,ij
2402  REAL :: x(nb(1),nb(2),nb(3)),y(nbo)
2403!---------------------------------------------------------------------
2404  ij = 0
2405  DO k=1,nb(3)
2406    DO j=1,nb(2)
2407      DO i=1,nb(1)
2408        ij = ij+1
2409        y(ij) = ATAN(x(i,j,k))
2410      ENDDO
2411    ENDDO
2412  ENDDO
2413!-
2414  nbo = nb(1)*nb(2)*nb(3)
2415  ma_atan_r31 = 0
2416!-----------------------
2417  END FUNCTION ma_atan_r31
2418!===
2419INTEGER FUNCTION ma_exp_r31(nb,x,nbo,y)
2420!---------------------------------------------------------------------
2421  IMPLICIT NONE
2422!-
2423  INTEGER :: nb(3),nbo,i,j,k,ij
2424  REAL :: x(nb(1),nb(2),nb(3)),y(nbo)
2425!---------------------------------------------------------------------
2426  ij = 0
2427  DO k=1,nb(3)
2428    DO j=1,nb(2)
2429      DO i=1,nb(1)
2430        ij = ij+1
2431        y(ij) = EXP(x(i,j,k))
2432      ENDDO
2433    ENDDO
2434  ENDDO
2435!-
2436  nbo = nb(1)*nb(2)*nb(3)
2437  ma_exp_r31 = 0
2438!----------------------
2439END FUNCTION ma_exp_r31
2440!===
2441INTEGER FUNCTION ma_log_r31(nb,x,nbo,y)
2442!---------------------------------------------------------------------
2443  IMPLICIT NONE
2444!-
2445  INTEGER :: nb(3),nbo,i,j,k,ij
2446  REAL :: x(nb(1),nb(2),nb(3)),y(nbo)
2447!---------------------------------------------------------------------
2448  ij = 0
2449  DO k=1,nb(3)
2450    DO j=1,nb(2)
2451      DO i=1,nb(1)
2452        ij = ij+1
2453        y(ij) = LOG(x(i,j,k))
2454      ENDDO
2455    ENDDO
2456  ENDDO
2457!-
2458  nbo = nb(1)*nb(2)*nb(3)
2459  ma_log_r31 = 0
2460!----------------------
2461END FUNCTION ma_log_r31
2462!===
2463INTEGER FUNCTION ma_sqrt_r31(nb,x,nbo,y)
2464!---------------------------------------------------------------------
2465  IMPLICIT NONE
2466!-
2467  INTEGER :: nb(3),nbo,i,j,k,ij
2468  REAL :: x(nb(1),nb(2),nb(3)),y(nbo)
2469!---------------------------------------------------------------------
2470  ij = 0
2471  DO k=1,nb(3)
2472    DO j=1,nb(2)
2473      DO i=1,nb(1)
2474        ij = ij+1
2475        y(ij) = SQRT(x(i,j,k))
2476      ENDDO
2477    ENDDO
2478  ENDDO
2479!-
2480  nbo = nb(1)*nb(2)*nb(3)
2481  ma_sqrt_r31 = 0
2482!-----------------------
2483END FUNCTION ma_sqrt_r31
2484!===
2485INTEGER FUNCTION ma_abs_r31(nb,x,nbo,y)
2486!---------------------------------------------------------------------
2487  IMPLICIT NONE
2488!-
2489  INTEGER :: nb(3),nbo,i,j,k,ij
2490  REAL :: x(nb(1),nb(2),nb(3)),y(nbo)
2491!---------------------------------------------------------------------
2492  ij = 0
2493  DO k=1,nb(3)
2494    DO j=1,nb(2)
2495      DO i=1,nb(1)
2496        ij = ij+1
2497        y(ij) = ABS(x(i,j,k))
2498      ENDDO
2499    ENDDO
2500  ENDDO
2501!-
2502  nbo = nb(1)*nb(2)*nb(3)
2503  ma_abs_r31 = 0
2504!----------------------
2505END FUNCTION ma_abs_r31
2506!===
2507INTEGER FUNCTION ma_chs_r31(nb,x,nbo,y)
2508!---------------------------------------------------------------------
2509  IMPLICIT NONE
2510!-
2511  INTEGER :: nb(3),nbo,i,j,k,ij
2512  REAL :: x(nb(1),nb(2),nb(3)),y(nbo)
2513!---------------------------------------------------------------------
2514  ij = 0
2515  DO k=1,nb(3)
2516    DO j=1,nb(2)
2517      DO i=1,nb(1)
2518        ij = ij+1
2519        y(ij) = x(i,j,k)*(-1.)
2520      ENDDO
2521    ENDDO
2522  ENDDO
2523!-
2524  nbo = nb(1)*nb(2)*nb(3)
2525  ma_chs_r31 = 0
2526!----------------------
2527END FUNCTION ma_chs_r31
2528!===
2529INTEGER FUNCTION ma_cels_r31(nb,x,nbo,y)
2530!---------------------------------------------------------------------
2531  IMPLICIT NONE
2532!-
2533  INTEGER :: nb(3),nbo,i,j,k,ij
2534  REAL :: x(nb(1),nb(2),nb(3)),y(nbo)
2535!---------------------------------------------------------------------
2536  ij = 0
2537  DO k=1,nb(3)
2538    DO j=1,nb(2)
2539      DO i=1,nb(1)
2540        ij = ij+1
2541        y(ij) = x(i,j,k)-273.15
2542      ENDDO
2543    ENDDO
2544  ENDDO
2545!-
2546  nbo = nb(1)*nb(2)*nb(3)
2547  ma_cels_r31 = 0
2548!-----------------------
2549END FUNCTION ma_cels_r31
2550!===
2551INTEGER FUNCTION ma_kelv_r31(nb,x,nbo,y)
2552!---------------------------------------------------------------------
2553  IMPLICIT NONE
2554!-
2555  INTEGER :: nb(3),nbo,i,j,k,ij
2556  REAL :: x(nb(1),nb(2),nb(3)),y(nbo)
2557!---------------------------------------------------------------------
2558  ij = 0
2559  DO k=1,nb(3)
2560    DO j=1,nb(2)
2561      DO i=1,nb(1)
2562        ij = ij+1
2563        y(ij) = x(i,j,k)+273.15
2564      ENDDO
2565    ENDDO
2566  ENDDO
2567!-
2568  nbo = nb(1)*nb(2)*nb(3)
2569  ma_kelv_r31 = 0
2570!-----------------------
2571  END FUNCTION ma_kelv_r31
2572!===
2573INTEGER FUNCTION ma_deg_r31(nb,x,nbo,y)
2574!---------------------------------------------------------------------
2575  IMPLICIT NONE
2576!-
2577  INTEGER :: nb(3),nbo,i,j,k,ij
2578  REAL :: x(nb(1),nb(2),nb(3)),y(nbo)
2579!---------------------------------------------------------------------
2580  ij = 0
2581  DO k=1,nb(3)
2582    DO j=1,nb(2)
2583      DO i=1,nb(1)
2584        ij = ij+1
2585        y(ij) = x(i,j,k)*57.29577951
2586      ENDDO
2587    ENDDO
2588  ENDDO
2589!-
2590  nbo = nb(1)*nb(2)*nb(3)
2591  ma_deg_r31 = 0
2592!----------------------
2593END FUNCTION ma_deg_r31
2594!===
2595INTEGER FUNCTION ma_rad_r31(nb,x,nbo,y)
2596!---------------------------------------------------------------------
2597  IMPLICIT NONE
2598!-
2599  INTEGER :: nb(3),nbo,i,j,k,ij
2600  REAL :: x(nb(1),nb(2),nb(3)),y(nbo)
2601!---------------------------------------------------------------------
2602  ij = 0
2603  DO k=1,nb(3)
2604    DO j=1,nb(2)
2605      DO i=1,nb(1)
2606        ij = ij+1
2607        y(ij) = x(i,j,k)*0.01745329252
2608      ENDDO
2609    ENDDO
2610  ENDDO
2611!-
2612  nbo = nb(1)*nb(2)*nb(3)
2613  ma_rad_r31 = 0
2614!----------------------
2615END FUNCTION ma_rad_r31
2616!===
2617INTEGER FUNCTION ma_ident_r31(nb,x,nbo,y)
2618!---------------------------------------------------------------------
2619  IMPLICIT NONE
2620!-
2621  INTEGER :: nb(3),nbo,i,j,k,ij
2622  REAL :: x(nb(1),nb(2),nb(3)),y(nbo)
2623!---------------------------------------------------------------------
2624  ij = 0
2625  DO k=1,nb(3)
2626    DO j=1,nb(2)
2627      DO i=1,nb(1)
2628        ij = ij+1
2629        y(ij) = x(i,j,k)
2630      ENDDO
2631    ENDDO
2632  ENDDO
2633!-
2634  nbo = nb(1)*nb(2)*nb(3)
2635  ma_ident_r31 = 0
2636!------------------------
2637END FUNCTION ma_ident_r31
2638!-
2639!=== OPERATIONS (two argument)
2640!-
2641INTEGER FUNCTION ma_add_r31(nb,x,s,nbo,y)
2642!---------------------------------------------------------------------
2643  IMPLICIT NONE
2644!-
2645  INTEGER :: nb(3),nbo
2646  REAL :: x(nb(1),nb(2),nb(3)),s,y(nbo)
2647!-
2648  INTEGER :: i,j,k,ij
2649!---------------------------------------------------------------------
2650  ij = 0
2651  DO k=1,nb(3)
2652    DO j=1,nb(2)
2653      DO i=1,nb(1)
2654        ij = ij+1
2655        y(ij) = x(i,j,k)+s
2656      ENDDO
2657    ENDDO
2658  ENDDO
2659!-
2660  nbo = nb(1)*nb(2)*nb(3)
2661  ma_add_r31 = 0
2662!----------------------
2663END FUNCTION ma_add_r31
2664!===
2665INTEGER FUNCTION ma_sub_r31(nb,x,s,nbo,y)
2666!---------------------------------------------------------------------
2667  IMPLICIT NONE
2668!-
2669  INTEGER :: nb(3),nbo
2670  REAL :: x(nb(1),nb(2),nb(3)),s,y(nbo)
2671!-
2672  INTEGER :: i,j,k,ij
2673!---------------------------------------------------------------------
2674  ij = 0
2675  DO k=1,nb(3)
2676    DO j=1,nb(2)
2677      DO i=1,nb(1)
2678        ij = ij+1
2679        y(ij) = x(i,j,k)-s
2680      ENDDO
2681    ENDDO
2682  ENDDO
2683!-
2684  nbo = nb(1)*nb(2)*nb(3)
2685  ma_sub_r31 = 0
2686!----------------------
2687END FUNCTION ma_sub_r31
2688!===
2689INTEGER FUNCTION ma_subi_r31(nb,x,s,nbo,y)
2690!---------------------------------------------------------------------
2691  IMPLICIT NONE
2692!-
2693  INTEGER :: nb(3),nbo
2694  REAL :: x(nb(1),nb(2),nb(3)),s,y(nbo)
2695!-
2696  INTEGER :: i,j,k,ij
2697!---------------------------------------------------------------------
2698  ij = 0
2699  DO k=1,nb(3)
2700    DO j=1,nb(2)
2701      DO i=1,nb(1)
2702        ij = ij+1
2703        y(ij) =  s-x(i,j,k)
2704      ENDDO
2705    ENDDO
2706  ENDDO
2707!-
2708  nbo = nb(1)*nb(2)*nb(3)
2709  ma_subi_r31 = 0
2710!-----------------------
2711END FUNCTION ma_subi_r31
2712!===
2713INTEGER FUNCTION ma_mult_r31(nb,x,s,nbo,y)
2714!---------------------------------------------------------------------
2715  IMPLICIT NONE
2716!-
2717  INTEGER :: nb(3),nbo
2718  REAL :: x(nb(1),nb(2),nb(3)),s,y(nbo)
2719!-
2720  INTEGER :: i,j,k,ij
2721!---------------------------------------------------------------------
2722  ij = 0
2723  DO k=1,nb(3)
2724    DO j=1,nb(2)
2725      DO i=1,nb(1)
2726        ij = ij+1
2727        y(ij) = x(i,j,k)*s
2728      ENDDO
2729    ENDDO
2730  ENDDO
2731!-
2732  nbo = nb(1)*nb(2)*nb(3)
2733  ma_mult_r31 = 0
2734!-----------------------
2735END FUNCTION ma_mult_r31
2736!===
2737INTEGER FUNCTION ma_div_r31(nb,x,s,nbo,y)
2738!---------------------------------------------------------------------
2739  IMPLICIT NONE
2740!-
2741  INTEGER :: nb(3),nbo
2742  REAL :: x(nb(1),nb(2),nb(3)),s,y(nbo)
2743!-
2744  INTEGER :: i,j,k,ij
2745!---------------------------------------------------------------------
2746  ij = 0
2747  DO k=1,nb(3)
2748    DO j=1,nb(2)
2749      DO i=1,nb(1)
2750        ij = ij+1
2751        y(ij) = x(i,j,k)/s
2752      ENDDO
2753    ENDDO
2754  ENDDO
2755!-
2756  nbo = nb(1)*nb(2)*nb(3)
2757  ma_div_r31 = 0
2758!----------------------
2759END FUNCTION ma_div_r31
2760!===
2761INTEGER FUNCTION ma_divi_r31(nb,x,s,nbo,y)
2762!---------------------------------------------------------------------
2763  IMPLICIT NONE
2764!-
2765  INTEGER :: nb(3),nbo
2766  REAL :: x(nb(1),nb(2),nb(3)),s,y(nbo)
2767!-
2768  INTEGER :: i,j,k,ij
2769!---------------------------------------------------------------------
2770  ij = 0
2771  DO k=1,nb(3)
2772    DO j=1,nb(2)
2773      DO i=1,nb(1)
2774        ij = ij+1
2775        y(ij) = s/x(i,j,k)
2776      ENDDO
2777    ENDDO
2778  ENDDO
2779!-
2780  nbo = nb(1)*nb(2)*nb(3)
2781  ma_divi_r31 = 0
2782!-----------------------
2783END FUNCTION ma_divi_r31
2784!===
2785INTEGER FUNCTION ma_power_r31(nb,x,s,nbo,y)
2786!---------------------------------------------------------------------
2787  IMPLICIT NONE
2788!-
2789  INTEGER :: nb(3),nbo
2790  REAL :: x(nb(1),nb(2),nb(3)),s,y(nbo)
2791!-
2792  INTEGER :: i,j,k,ij
2793!---------------------------------------------------------------------
2794  ij = 0
2795  DO k=1,nb(3)
2796    DO j=1,nb(2)
2797      DO i=1,nb(1)
2798        ij = ij+1
2799        y(ij) = x(i,j,k)**s
2800      ENDDO
2801    ENDDO
2802  ENDDO
2803!-
2804  nbo = nb(1)*nb(2)*nb(3)
2805  ma_power_r31 = 0
2806!------------------------
2807END FUNCTION ma_power_r31
2808!===
2809INTEGER FUNCTION ma_fumin_r31(nb,x,s,nbo,y)
2810!---------------------------------------------------------------------
2811  IMPLICIT NONE
2812!-
2813  INTEGER :: nb(3),nbo
2814  REAL :: x(nb(1),nb(2),nb(3)),s,y(nbo)
2815!-
2816  INTEGER :: i,j,k,ij
2817!---------------------------------------------------------------------
2818  ij = 0
2819  DO k=1,nb(3)
2820    DO j=1,nb(2)
2821      DO i=1,nb(1)
2822        ij = ij+1
2823        y(ij) = MIN(x(i,j,k),s)
2824      ENDDO
2825    ENDDO
2826  ENDDO
2827!-
2828  nbo = nb(1)*nb(2)*nb(3)
2829  ma_fumin_r31 = 0
2830!------------------------
2831END FUNCTION ma_fumin_r31
2832!===
2833INTEGER FUNCTION ma_fumax_r31(nb,x,s,nbo,y)
2834!---------------------------------------------------------------------
2835  IMPLICIT NONE
2836!-
2837  INTEGER :: nb(3),nbo
2838  REAL :: x(nb(1),nb(2),nb(3)),s,y(nbo)
2839!-
2840  INTEGER :: i,j,k,ij
2841!---------------------------------------------------------------------
2842  ij = 0
2843  DO k=1,nb(3)
2844    DO j=1,nb(2)
2845      DO i=1,nb(1)
2846        ij = ij+1
2847        y(ij) = MAX(x(i,j,k),s)
2848      ENDDO
2849    ENDDO
2850  ENDDO
2851!-
2852  nbo = nb(1)*nb(2)*nb(3)
2853  ma_fumax_r31 = 0
2854!------------------------
2855END FUNCTION ma_fumax_r31
2856!===
2857INTEGER FUNCTION ma_fuscat_r31(nb,x,nbi,ind,miss_val,nbo,y)
2858!---------------------------------------------------------------------
2859  IMPLICIT NONE
2860!-
2861  INTEGER :: nb(3),nbo,nbi
2862  INTEGER :: ind(nbi)
2863  REAL :: x(nb(1),nb(2),nb(3)),miss_val,y(nbo)
2864!-
2865  INTEGER :: i,j,k,ij,ii,ipos,ipp,isb
2866!---------------------------------------------------------------------
2867  ma_fuscat_r31 = 0
2868!-
2869  y(1:nbo) = miss_val
2870!-
2871  IF (nbi <= nb(1)*nb(2)*nb(3)) THEN
2872    ipos = 0
2873    isb = nb(1)*nb(2)
2874    DO ij=1,nbi
2875      IF (ind(ij) <= nbo .AND. ind(ij) > 0) THEN
2876        ipos = ipos+1
2877        k = ((ipos-1)/isb)+1
2878        ipp = ipos-(k-1)*isb
2879        j = ((ipp-1)/nb(1))+1
2880        i = (ipp-(j-1)*nb(1))
2881        y(ind(ij)) = x(i,j,k)
2882      ELSE
2883        IF (ind(ij) > nbo) ma_fuscat_r31  = ma_fuscat_r31+1
2884      ENDIF
2885    ENDDO
2886!-- Repeat the data if needed
2887    IF (MINVAL(ind) < 0) THEN
2888      DO i=1,nbi
2889        IF (ind(i) <= 0) THEN
2890          DO ii=1,ABS(ind(i))-1
2891            IF (ind(i+1)+ii <= nbo) THEN
2892              y(ind(i+1)+ii) = y(ind(i+1))
2893            ELSE
2894              ma_fuscat_r31  = ma_fuscat_r31+1
2895            ENDIF
2896          ENDDO
2897        ENDIF
2898      ENDDO
2899    ENDIF
2900  ELSE
2901    ma_fuscat_r31  = 1
2902  ENDIF
2903!-------------------------
2904END FUNCTION ma_fuscat_r31
2905!===
2906INTEGER FUNCTION ma_fugath_r31(nb,x,nbi,ind,miss_val,nbo,y)
2907!---------------------------------------------------------------------
2908  IMPLICIT NONE
2909!-
2910  INTEGER :: nb(3),nbo,nbi
2911  INTEGER :: ind(nbi)
2912  REAL :: x(nb(1),nb(2),nb(3)),miss_val,y(nbo)
2913!-
2914  INTEGER :: i,j,k,ij,ipos,ipp,isb
2915!---------------------------------------------------------------------
2916  IF (nbi <= nbo) THEN
2917    ma_fugath_r31 = 0
2918    y(1:nbo) = miss_val
2919    ipos = 0
2920    isb = nb(1)*nb(2)
2921    DO ij=1,nbi
2922      IF (ipos+1 <= nbo) THEN
2923        IF (ind(ij) > 0) THEN
2924          k = ((ind(ij)-1)/isb)+1
2925          ipp = ind(ij)-(k-1)*isb
2926          j = ((ipp-1)/nb(1))+1
2927          i = (ipp-(j-1)*nb(1))
2928          ipos = ipos+1
2929          y(ipos) = x(i,j,k)
2930        ENDIF
2931      ELSE
2932        IF (ipos+1 > nbo) ma_fugath_r31  = ma_fugath_r31+1
2933      ENDIF
2934    ENDDO
2935  ELSE
2936    ma_fugath_r31 = 1
2937  ENDIF
2938  nbo = ipos
2939!-------------------------
2940END FUNCTION ma_fugath_r31
2941!===
2942INTEGER FUNCTION ma_fufill_r31(nb,x,nbi,ind,miss_val,nbo,y)
2943!---------------------------------------------------------------------
2944  IMPLICIT NONE
2945!-
2946  INTEGER :: nb(3),nbo,nbi
2947  INTEGER :: ind(nbi)
2948  REAL :: x(nb(1),nb(2),nb(3)),miss_val,y(nbo)
2949!-
2950  INTEGER :: i,j,k,ij,ii,ipos,ipp,isb
2951!---------------------------------------------------------------------
2952  ma_fufill_r31 = 0
2953  IF (nbi <= nb(1)*nb(2)*nb(3)) THEN
2954    ipos = 0
2955    isb = nb(1)*nb(2)
2956    DO ij=1,nbi
2957      IF (ind(ij) <= nbo .AND. ind(ij) > 0) THEN
2958        ipos = ipos+1
2959        k = ((ipos-1)/isb)+1
2960        ipp = ipos-(k-1)*isb
2961        j = ((ipp-1)/nb(1))+1
2962        i = (ipp-(j-1)*nb(1))
2963        y(ind(ij)) = x(i,j,k)
2964      ELSE
2965        IF (ind(ij) > nbo) ma_fufill_r31  = ma_fufill_r31+1
2966      ENDIF
2967    ENDDO
2968!-- Repeat the data if needed
2969    IF (MINVAL(ind) < 0) THEN
2970      DO i=1,nbi
2971        IF (ind(i) <= 0) THEN
2972          DO ii=1,ABS(ind(i))-1
2973            IF (ind(i+1)+ii <= nbo) THEN
2974              y(ind(i+1)+ii) = y(ind(i+1))
2975            ELSE
2976              ma_fufill_r31  = ma_fufill_r31+1
2977            ENDIF
2978          ENDDO
2979        ENDIF
2980      ENDDO
2981    ENDIF
2982  ELSE
2983    ma_fufill_r31  = 1
2984  ENDIF
2985!-------------------------
2986END FUNCTION ma_fufill_r31
2987!===
2988INTEGER FUNCTION ma_fucoll_r31(nb,x,nbi,ind,miss_val,nbo,y)
2989!---------------------------------------------------------------------
2990  IMPLICIT NONE
2991!-
2992  INTEGER :: nb(3),nbo,nbi
2993  INTEGER :: ind(nbi)
2994  REAL :: x(nb(1),nb(2),nb(3)),miss_val,y(nbo)
2995!-
2996  INTEGER :: i,j,k,ij,ipos,ipp,isb
2997!---------------------------------------------------------------------
2998  IF (nbi <= nbo) THEN
2999    ma_fucoll_r31 = 0
3000    ipos = 0
3001    isb = nb(1)*nb(2)
3002    DO ij=1,nbi
3003      IF (ipos+1 <= nbo) THEN
3004        IF (ind(ij) > 0) THEN
3005          k = ((ind(ij)-1)/isb)+1
3006          ipp = ind(ij)-(k-1)*isb
3007          j = ((ipp-1)/nb(1))+1
3008          i = (ipp-(j-1)*nb(1))
3009          ipos = ipos+1
3010          y(ipos) = x(i,j,k)
3011        ENDIF
3012      ELSE
3013        IF (ipos+1 > nbo) ma_fucoll_r31  = ma_fucoll_r31+1
3014      ENDIF
3015    ENDDO
3016  ELSE
3017    ma_fucoll_r31 = 1
3018  ENDIF
3019  nbo = ipos
3020!-------------------------
3021END FUNCTION ma_fucoll_r31
3022!===
3023INTEGER FUNCTION ma_fuundef_r31(nb,x,nbi,ind,miss_val,nbo,y)
3024!---------------------------------------------------------------------
3025  IMPLICIT NONE
3026!-
3027  INTEGER :: nb(3),nbo,nbi
3028  INTEGER :: ind(nbi)
3029  REAL :: x(nb(1),nb(2),nb(3)),miss_val,y(nbo)
3030!-
3031  INTEGER :: i,j,k,ij,ipp,isb
3032!---------------------------------------------------------------------
3033  IF (nbi <= nbo .AND. nbo == nb(1)*nb(2)*nb(3)) THEN
3034    ma_fuundef_r31 = 0
3035    isb = nb(1)*nb(2)
3036    DO ij=1,nbo
3037      k = ((ij-1)/isb)+1
3038      ipp = ij-(k-1)*isb
3039      j = ((ipp-1)/nb(1))+1
3040      i = (ipp-(j-1)*nb(1))
3041      y(ij) = x(i,j,k)
3042    ENDDO
3043    DO i=1,nbi
3044      IF (ind(i) <= nbo .AND. ind(i) > 0) THEN
3045        y(ind(i)) =  miss_val
3046      ELSE
3047        IF (ind(i) > nbo) ma_fuundef_r31  = ma_fuundef_r31+1
3048      ENDIF
3049    ENDDO
3050  ELSE
3051    ma_fuundef_r31 = 1
3052  ENDIF
3053!--------------------------
3054END FUNCTION ma_fuundef_r31
3055!===
3056INTEGER FUNCTION ma_fuonly_r31(nb,x,nbi,ind,miss_val,nbo,y)
3057!---------------------------------------------------------------------
3058  IMPLICIT NONE
3059!-
3060  INTEGER :: nb(3),nbo,nbi
3061  INTEGER :: ind(nbi)
3062  REAL :: x(nb(1),nb(2),nb(3)),miss_val,y(nbo)
3063!-
3064  INTEGER :: i,j,k,ij,ipp,isb
3065!---------------------------------------------------------------------
3066  IF (     (nbi <= nbo).AND.(nbo == nb(1)*nb(2)*nb(3)) &
3067 &    .AND.ALL(ind(1:nbi) <= nbo) ) THEN
3068    ma_fuonly_r31 = 0
3069    y(1:nbo) = miss_val
3070    isb = nb(1)*nb(2)
3071    DO ij=1,nbi
3072      IF (ind(ij) > 0) THEN
3073        k = ((ind(ij)-1)/isb)+1
3074        ipp = ind(ij)-(k-1)*isb
3075        j = ((ipp-1)/nb(1))+1
3076        i = (ipp-(j-1)*nb(1))
3077        y(ind(ij)) = x(i,j,k)
3078      ENDIF
3079    ENDDO
3080  ELSE
3081    ma_fuonly_r31 = 1
3082  ENDIF
3083!-------------------------
3084END FUNCTION ma_fuonly_r31
3085!===
3086SUBROUTINE moycum (opp,np,px,py,pwx)
3087!---------------------------------------------------------------------
3088!- Does time operations
3089!---------------------------------------------------------------------
3090  IMPLICIT NONE
3091!-
3092  CHARACTER(LEN=7) :: opp
3093  INTEGER :: np
3094  REAL,DIMENSION(:) :: px,py
3095  INTEGER :: pwx
3096!---------------------------------------------------------------------
3097  IF (pwx /= 0) THEN
3098    IF      (opp == 'ave') THEN
3099      px(1:np)=(px(1:np)*pwx+py(1:np))/REAL(pwx+1)
3100    ELSE IF (opp == 't_sum') THEN
3101      px(1:np)=px(1:np)+py(1:np)
3102    ELSE IF ( (opp == 'l_min').OR.(opp == 't_min') ) THEN
3103      px(1:np)=MIN(px(1:np),py(1:np))
3104    ELSE IF ( (opp == 'l_max').OR.(opp == 't_max') ) THEN
3105      px(1:np)=MAX(px(1:np),py(1:np))
3106    ELSE
3107      CALL ipslerr(3,"moycum",'Unknown time operation',opp,' ')
3108    ENDIF
3109  ELSE
3110    IF      (opp == 'l_min') THEN
3111      px(1:np)=MIN(px(1:np),py(1:np))
3112    ELSE IF (opp == 'l_max') THEN
3113      px(1:np)=MAX(px(1:np),py(1:np))
3114    ELSE
3115      px(1:np)=py(1:np)
3116    ENDIF
3117  ENDIF
3118!--------------------
3119END SUBROUTINE moycum
3120!===
3121SUBROUTINE moycum_index( opp, px, py, pwx, nbi, ind )
3122!---------------------------------------------------------------------
3123!- Does time operations on index points
3124!---------------------------------------------------------------------
3125  IMPLICIT NONE
3126!-
3127
3128  !! 0. Parameters and variables declaration
3129
3130  !! 0.1 Input variables
3131
3132  CHARACTER(LEN=7), INTENT(in)        :: opp  !! Operation performed
3133  INTEGER, INTENT(in)                 :: nbi  !! Size of index vector
3134  INTEGER, DIMENSION(nbi), INTENT(in) :: ind  !! Index vector
3135  REAL, DIMENSION(:), INTENT(in)      :: py   !! Vector containing the
3136                                              !! previous values of px
3137                                              !! Warning : due to memory
3138                                              !! optimization, we have
3139                                              !! generally SIZE(px) /= SIZE(py) 
3140  INTEGER, INTENT(in)                 :: pwx  !! Used to calculate average value                             
3141
3142  !! 0.3 Modified variables
3143
3144  REAL, DIMENSION(:), INTENT(inout)   :: px   !! Result
3145
3146  !! 0.4 Local variables
3147
3148  INTEGER :: ig                               !! Index
3149
3150!---------------------------------------------------------------------
3151
3152  !! Perform operations only if the values of ind don't exceed the size of px
3153
3154  IF ( MAXVAL(ind) > SIZE(px) ) THEN
3155     CALL ipslerr(3,"moycum_index", &
3156          & "the index vector is out of range for px", &
3157          & "Indexation vector problem. We stop", " " )
3158  END IF
3159
3160  IF (pwx /= 0) THEN
3161     IF      (opp == 'ave') THEN
3162        DO ig = 1,nbi
3163           px(ind(ig)) = (px(ind(ig))*pwx + py(ind(ig)))/REAL(pwx+1)
3164        END DO
3165     ELSE IF (opp == 't_sum') THEN
3166        DO ig = 1,nbi
3167           px(ind(ig)) = px(ind(ig)) + py(ind(ig))
3168        END DO
3169     ELSE IF ( (opp == 'l_min').OR.(opp == 't_min') ) THEN
3170        DO ig = 1,nbi
3171           px(ind(ig)) = MIN(px(ind(ig)),py(ind(ig)))
3172        END DO
3173     ELSE IF ( (opp == 'l_max').OR.(opp == 't_max') ) THEN
3174        DO ig = 1,nbi
3175           px(ind(ig)) = MAX(px(ind(ig)),py(ind(ig)))
3176        END DO
3177     ELSE
3178        CALL ipslerr(3,"moycum_index",'Unknown time operation',opp,' ')
3179     END IF
3180  ELSE
3181    IF      (opp == 'l_min') THEN
3182       DO ig = 1,nbi
3183          px(ind(ig)) = MIN(px(ind(ig)),py(ind(ig)))
3184       END DO
3185    ELSE IF (opp == 'l_max') THEN
3186       DO ig = 1,nbi
3187          px(ind(ig)) = MAX(px(ind(ig)),py(ind(ig)))
3188       END DO
3189    ELSE
3190       DO ig = 1,nbi
3191          px(ind(ig)) = py(ind(ig))
3192       END DO
3193    ENDIF
3194 END IF
3195
3196END SUBROUTINE moycum_index
3197
3198
3199!-----------------
3200END MODULE mathelp
Note: See TracBrowser for help on using the repository browser.