source: vendors/IOIPSL/current/src/mathelp.f90 @ 1895

Last change on this file since 1895 was 1895, checked in by flavoni, 11 years ago

importing IOIPSL on vendors

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