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

Last change on this file since 11 was 11, checked in by bellier, 17 years ago

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