source: branches/UKMO/r6232_tracer_advection/NEMOGCM/EXTERNAL/IOIPSL/src/mathelp.f90 @ 9295

Last change on this file since 9295 was 9295, checked in by jcastill, 3 years ago

Remove svn keywords

File size: 80.4 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,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!===
3121!-----------------
3122END MODULE mathelp
Note: See TracBrowser for help on using the repository browser.