New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
mathelp.f90 in utils/tools_dev_r12970_AGRIF_CMEMS/DOMAINcfg/src – NEMO

source: utils/tools_dev_r12970_AGRIF_CMEMS/DOMAINcfg/src/mathelp.f90 @ 13024

Last change on this file since 13024 was 13024, checked in by rblod, 4 years ago

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