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

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 30 - (show annotations)
Thu Apr 1 09:07:28 2010 UTC (14 years, 2 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 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