/[lmdze]/trunk/libf/IOIPSL/mathelp.f90
ViewVC logotype

Annotation of /trunk/libf/IOIPSL/mathelp.f90

Parent Directory Parent Directory | Revision Log Revision Log

Revision 30 - (hide annotations)
Thu Apr 1 09:07:28 2010 UTC (14 years, 3 months ago) by guez
File size: 88787 byte(s) 
Imported Source files of the external library "IOIPSL_Lionel" into
"libf/IOIPSL".

Split "cray.f90" into "scopy.f90" and "ssum.f90".

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