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

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 32 - (hide annotations)
Tue Apr 6 17:52:58 2010 UTC (14 years, 1 month ago) by guez
File size: 71074 byte(s) 
Split "stringop.f90" into single-procedure files. Gathered files in directory
"IOIPSL/Stringop".

Split "flincom.f90" into "flincom.f90" and "flinget.f90". Removed
unused procedures from module "flincom". Removed unused argument
"filename" of procedure "flinopen_nozoom".

Removed unused files.

Split "grid_change.f90" into "grid_change.f90" and
"gr_phy_write_3d.f90".

Removed unused procedures from modules "calendar", "ioipslmpp",
"grid_atob", "gath_cpl" and "getincom". Removed unused procedures in
files "ppm3d.f" and "thermcell.f".

Split "mathelp.f90" into "mathelp.f90" and "mathop.f90".

Removed unused variable "dpres" of module "comvert".

Use argument "itau" instead of local variables "iadvtr" and "first" to
control algorithm in procedure "fluxstokenc".

Removed unused arguments of procedure "integrd".

Removed useless computations at the end of procedure "leapfrog".

Merged common block "matrfil" into module "parafilt".
1 guez 32 MODULE mathop_m
2    
3     ! From mathelp.f90, version 2.0 2004/04/05 14:47:50
4    
5     implicit none
6    
7     PRIVATE
8     PUBLIC :: mathop
9    
10     INTERFACE mathop
11     MODULE PROCEDURE mathop_r11, mathop_r21, mathop_r31
12     END INTERFACE
13    
14     !- Variables used to detect and identify the operations
15     CHARACTER(LEN=120), SAVE :: indexfu = &
16     'gather, scatter, fill, coll, undef, only'
17    
18     CONTAINS
19    
20     SUBROUTINE mathop_r11 &
21     & (fun, nb, work_in, miss_val, nb_index, nindex, scal, nb_max, work_out)
22     !- This subroutines gives an interface to the various operation
23     !- which are allowed. The interface is general enough to allow its use
24     !- for other cases.
25    
26     !- INPUT
27    
28     !- fun : function to be applied to the vector of data
29     !- nb : Length of input vector
30     !- work_in : Input vector of data (REAL)
31     !- miss_val : The value of the missing data flag (it has to be a
32     !- maximum value, in f90 : huge( a real ))
33     !- nb_index : Length of index vector
34     !- nindex : Vector of indices
35     !- scal : A scalar value for vector/scalar operations
36     !- nb_max : maximum length of output vector
37    
38     !- OUTPUT
39    
40     !- nb_max : Actual length of output variable
41     !- work_out : Output vector after the operation was applied
42     USE errioipsl, ONLY : histerr
43    
44     CHARACTER(LEN=7) :: fun
45     INTEGER :: nb, nb_max, nb_index
46     INTEGER :: nindex(nb_index)
47     REAL :: work_in(nb), scal, miss_val
48     REAL :: work_out(nb_max)
49    
50     INTEGER :: ierr
51    
52     INTRINSIC SIN, COS, TAN, ASIN, ACOS, ATAN, EXP, ALOG, SQRT, ABS
53     !---------------------------------------------------------------------
54     ierr = 0
55    
56     IF (scal >= miss_val-1.) THEN
57     IF (INDEX(indexfu, fun(1:LEN_TRIM(fun))) == 0) THEN
58     SELECT CASE (fun)
59     CASE('sin')
60     ierr = ma_sin_r11(nb, work_in, nb_max, work_out)
61     CASE('cos')
62     ierr = ma_cos_r11(nb, work_in, nb_max, work_out)
63     CASE('tan')
64     ierr = ma_tan_r11(nb, work_in, nb_max, work_out)
65     CASE('asin')
66     ierr = ma_asin_r11(nb, work_in, nb_max, work_out)
67     CASE('acos')
68     ierr = ma_acos_r11(nb, work_in, nb_max, work_out)
69     CASE('atan')
70     ierr = ma_atan_r11(nb, work_in, nb_max, work_out)
71     CASE('exp')
72     ierr = ma_exp_r11(nb, work_in, nb_max, work_out)
73     CASE('log')
74     ierr = ma_alog_r11(nb, work_in, nb_max, work_out)
75     CASE('sqrt')
76     ierr = ma_sqrt_r11(nb, work_in, nb_max, work_out)
77     CASE('chs')
78     ierr = ma_chs_r11(nb, work_in, nb_max, work_out)
79     CASE('abs')
80     ierr = ma_abs_r11(nb, work_in, nb_max, work_out)
81     CASE('cels')
82     ierr = ma_cels_r11(nb, work_in, nb_max, work_out)
83     CASE('kelv')
84     ierr = ma_kelv_r11(nb, work_in, nb_max, work_out)
85     CASE('deg')
86     ierr = ma_deg_r11(nb, work_in, nb_max, work_out)
87     CASE('rad')
88     ierr = ma_rad_r11(nb, work_in, nb_max, work_out)
89     CASE('ident')
90     ierr = ma_ident_r11(nb, work_in, nb_max, work_out)
91     CASE DEFAULT
92     CALL histerr(3, "mathop", &
93     & 'scalar variable undefined and no indexing', &
94     & 'but still unknown function', fun)
95     END SELECT
96     IF (ierr > 0) THEN
97     CALL histerr(3, "mathop", &
98     & 'Error while executing a simple function', fun, ' ')
99     ENDIF
100     ELSE
101     SELECT CASE (fun)
102     CASE('gather')
103     ierr = ma_fugath_r11(nb, work_in, nb_index, nindex, &
104     & miss_val, nb_max, work_out)
105     CASE('scatter')
106     IF (nb_index > nb) THEN
107     work_out(1:nb_max) = miss_val
108     ierr=1
109     ELSE
110     ierr = ma_fuscat_r11(nb, work_in, nb_index, nindex, &
111     & miss_val, nb_max, work_out)
112     ENDIF
113     CASE('coll')
114     ierr = ma_fucoll_r11(nb, work_in, nb_index, nindex, &
115     & miss_val, nb_max, work_out)
116     CASE('fill')
117     ierr = ma_fufill_r11(nb, work_in, nb_index, nindex, &
118     & miss_val, nb_max, work_out)
119     CASE('undef')
120     ierr = ma_fuundef_r11(nb, work_in, nb_index, nindex, &
121     & miss_val, nb_max, work_out)
122     CASE('only')
123     ierr = ma_fuonly_r11(nb, work_in, nb_index, nindex, &
124     & miss_val, nb_max, work_out)
125     CASE DEFAULT
126     CALL histerr(3, "mathop", &
127     & 'scalar variable undefined and indexing', &
128     & 'was requested but with unknown function', fun)
129     END SELECT
130     IF (ierr > 0) THEN
131     CALL histerr(3, "mathop_r11", &
132     & 'Error while executing an indexing function', fun, ' ')
133     ENDIF
134     ENDIF
135     ELSE
136     SELECT CASE (fun)
137     CASE('fumin')
138     ierr = ma_fumin_r11(nb, work_in, scal, nb_max, work_out)
139     CASE('fumax')
140     ierr = ma_fumax_r11(nb, work_in, scal, nb_max, work_out)
141     CASE('add')
142     ierr = ma_add_r11(nb, work_in, scal, nb_max, work_out)
143     CASE('subi')
144     ierr = ma_subi_r11(nb, work_in, scal, nb_max, work_out)
145     CASE('sub')
146     ierr = ma_sub_r11(nb, work_in, scal, nb_max, work_out)
147     CASE('mult')
148     ierr = ma_mult_r11(nb, work_in, scal, nb_max, work_out)
149     CASE('div')
150     ierr = ma_div_r11(nb, work_in, scal, nb_max, work_out)
151     CASE('divi')
152     ierr = ma_divi_r11(nb, work_in, scal, nb_max, work_out)
153     CASE('power')
154     ierr = ma_power_r11(nb, work_in, scal, nb_max, work_out)
155     CASE DEFAULT
156     CALL histerr(3, "mathop", &
157     & 'Unknown operation with a scalar', fun, ' ')
158     END SELECT
159     IF (ierr > 0) THEN
160     CALL histerr(3, "mathop", &
161     & 'Error while executing a scalar function', fun, ' ')
162     ENDIF
163     ENDIF
164     !------------------------
165     END SUBROUTINE mathop_r11
166     !-
167     !=== FUNCTIONS (only one argument)
168     !-
169     INTEGER FUNCTION ma_sin_r11(nb, x, nbo, y)
170     INTEGER :: nb, nbo, i
171     REAL :: x(nb), y(nbo)
172     !---------------------------------------------------------------------
173     DO i=1, nb
174     y(i) = SIN(x(i))
175     ENDDO
176    
177     nbo = nb
178     ma_sin_r11 = 0
179     !----------------------
180     END FUNCTION ma_sin_r11
181    
182     !************************************************
183    
184     INTEGER FUNCTION ma_cos_r11(nb, x, nbo, y)
185     INTEGER :: nb, nbo, i
186     REAL :: x(nb), y(nbo)
187     !---------------------------------------------------------------------
188     DO i=1, nb
189     y(i) = COS(x(i))
190     ENDDO
191    
192     nbo = nb
193     ma_cos_r11 = 0
194     !----------------------
195     END FUNCTION ma_cos_r11
196    
197     !************************************************
198    
199     INTEGER FUNCTION ma_tan_r11(nb, x, nbo, y)
200     INTEGER :: nb, nbo, i
201     REAL :: x(nb), y(nbo)
202     !---------------------------------------------------------------------
203     DO i=1, nb
204     y(i) = TAN(x(i))
205     ENDDO
206    
207     nbo = nb
208     ma_tan_r11 = 0
209     !----------------------
210     END FUNCTION ma_tan_r11
211    
212     !************************************************
213    
214     INTEGER FUNCTION ma_asin_r11(nb, x, nbo, y)
215     INTEGER :: nb, nbo, i
216     REAL :: x(nb), y(nbo)
217     !---------------------------------------------------------------------
218     DO i=1, nb
219     y(i) = ASIN(x(i))
220     ENDDO
221    
222     nbo = nb
223     ma_asin_r11 = 0
224     !-----------------------
225     END FUNCTION ma_asin_r11
226    
227     !************************************************
228    
229     INTEGER FUNCTION ma_acos_r11(nb, x, nbo, y)
230     INTEGER :: nb, nbo, i
231     REAL :: x(nb), y(nbo)
232     !---------------------------------------------------------------------
233     DO i=1, nb
234     y(i) = ACOS(x(i))
235     ENDDO
236    
237     nbo = nb
238     ma_acos_r11 = 0
239     !-----------------------
240     END FUNCTION ma_acos_r11
241    
242     !************************************************
243    
244     INTEGER FUNCTION ma_atan_r11(nb, x, nbo, y)
245     INTEGER :: nb, nbo, i
246     REAL :: x(nb), y(nbo)
247     !---------------------------------------------------------------------
248     DO i=1, nb
249     y(i) = ATAN(x(i))
250     ENDDO
251    
252     nbo = nb
253     ma_atan_r11 = 0
254     !-----------------------
255     END FUNCTION ma_atan_r11
256    
257     !************************************************
258    
259     INTEGER FUNCTION ma_exp_r11(nb, x, nbo, y)
260     INTEGER :: nb, nbo, i
261     REAL :: x(nb), y(nbo)
262     !---------------------------------------------------------------------
263     DO i=1, nb
264     y(i) = EXP(x(i))
265     ENDDO
266    
267     nbo = nb
268     ma_exp_r11 = 0
269     !----------------------
270     END FUNCTION ma_exp_r11
271    
272     !************************************************
273    
274     INTEGER FUNCTION ma_alog_r11(nb, x, nbo, y)
275     INTEGER :: nb, nbo, i
276     REAL :: x(nb), y(nbo)
277     !---------------------------------------------------------------------
278     DO i=1, nb
279     y(i) = alog(x(i))
280     ENDDO
281    
282     nbo = nb
283     ma_alog_r11 = 0
284     !-----------------------
285     END FUNCTION ma_alog_r11
286    
287     !************************************************
288    
289     INTEGER FUNCTION ma_sqrt_r11(nb, x, nbo, y)
290     INTEGER :: nb, nbo, i
291     REAL :: x(nb), y(nbo)
292     !---------------------------------------------------------------------
293     DO i=1, nb
294     y(i) = SQRT(x(i))
295     ENDDO
296    
297     nbo = nb
298     ma_sqrt_r11 = 0
299     !-----------------------
300     END FUNCTION ma_sqrt_r11
301    
302     !************************************************
303    
304     INTEGER FUNCTION ma_abs_r11(nb, x, nbo, y)
305     INTEGER :: nb, nbo, i
306     REAL :: x(nb), y(nbo)
307     !---------------------------------------------------------------------
308     DO i=1, nb
309     y(i) = ABS(x(i))
310     ENDDO
311    
312     nbo = nb
313     ma_abs_r11 = 0
314     !----------------------
315     END FUNCTION ma_abs_r11
316    
317     !************************************************
318    
319     INTEGER FUNCTION ma_chs_r11(nb, x, nbo, y)
320     INTEGER :: nb, nbo, i
321     REAL :: x(nb), y(nbo)
322     !---------------------------------------------------------------------
323     DO i=1, nb
324     y(i) = x(i)*(-1.)
325     ENDDO
326    
327     nbo = nb
328     ma_chs_r11 = 0
329     !----------------------
330     END FUNCTION ma_chs_r11
331    
332     !************************************************
333    
334     INTEGER FUNCTION ma_cels_r11(nb, x, nbo, y)
335     INTEGER :: nb, nbo, i
336     REAL :: x(nb), y(nbo)
337     !---------------------------------------------------------------------
338     DO i=1, nb
339     y(i) = x(i)-273.15
340     ENDDO
341    
342     nbo = nb
343     ma_cels_r11 = 0
344     !-----------------------
345     END FUNCTION ma_cels_r11
346    
347     !************************************************
348    
349     INTEGER FUNCTION ma_kelv_r11(nb, x, nbo, y)
350     INTEGER :: nb, nbo, i
351     REAL :: x(nb), y(nbo)
352     !---------------------------------------------------------------------
353     DO i=1, nb
354     y(i) = x(i)+273.15
355     ENDDO
356    
357     nbo = nb
358     ma_kelv_r11 = 0
359     !-----------------------
360     END FUNCTION ma_kelv_r11
361    
362     !************************************************
363    
364     INTEGER FUNCTION ma_deg_r11(nb, x, nbo, y)
365     INTEGER :: nb, nbo, i
366     REAL :: x(nb), y(nbo)
367     !---------------------------------------------------------------------
368     DO i=1, nb
369     y(i) = x(i)*57.29577951
370     ENDDO
371    
372     nbo = nb
373     ma_deg_r11 = 0
374     !-----------------------
375     END FUNCTION ma_deg_r11
376    
377     !************************************************
378    
379     INTEGER FUNCTION ma_rad_r11(nb, x, nbo, y)
380     INTEGER :: nb, nbo, i
381     REAL :: x(nb), y(nbo)
382     !---------------------------------------------------------------------
383     DO i=1, nb
384     y(i) = x(i)*0.01745329252
385     ENDDO
386    
387     nbo = nb
388     ma_rad_r11 = 0
389     !----------------------
390     END FUNCTION ma_rad_r11
391    
392     !************************************************
393    
394     INTEGER FUNCTION ma_ident_r11(nb, x, nbo, y)
395     INTEGER :: nb, nbo, i
396     REAL :: x(nb), y(nbo)
397     !---------------------------------------------------------------------
398     DO i=1, nb
399     y(i) = x(i)
400     ENDDO
401    
402     nbo = nb
403     ma_ident_r11 = 0
404     !------------------------
405     END FUNCTION ma_ident_r11
406     !-
407     !=== OPERATIONS (two argument)
408     !-
409     INTEGER FUNCTION ma_add_r11(nb, x, s, nbo, y)
410     INTEGER :: nb, nbo
411     REAL :: x(nb), s, y(nbo)
412    
413     INTEGER :: i
414     !---------------------------------------------------------------------
415     DO i=1, nb
416     y(i) = x(i)+s
417     ENDDO
418    
419     nbo = nb
420     ma_add_r11 = 0
421     !-----------------------
422     END FUNCTION ma_add_r11
423    
424     !************************************************
425    
426     INTEGER FUNCTION ma_sub_r11(nb, x, s, nbo, y)
427     INTEGER :: nb, nbo
428     REAL :: x(nb), s, y(nbo)
429    
430     INTEGER :: i
431     !---------------------------------------------------------------------
432     DO i=1, nb
433     y(i) = x(i)-s
434     ENDDO
435    
436     nbo = nb
437     ma_sub_r11 = 0
438     !----------------------
439     END FUNCTION ma_sub_r11
440    
441     !************************************************
442    
443     INTEGER FUNCTION ma_subi_r11(nb, x, s, nbo, y)
444     INTEGER :: nb, nbo
445     REAL :: x(nb), s, y(nbo)
446    
447     INTEGER :: i
448     !---------------------------------------------------------------------
449     DO i=1, nb
450     y(i) = s-x(i)
451     ENDDO
452    
453     nbo = nb
454     ma_subi_r11 = 0
455     !-----------------------
456     END FUNCTION ma_subi_r11
457    
458     !************************************************
459    
460     INTEGER FUNCTION ma_mult_r11(nb, x, s, nbo, y)
461     INTEGER :: nb, nbo
462     REAL :: x(nb), s, y(nbo)
463    
464     INTEGER :: i
465     !---------------------------------------------------------------------
466     DO i=1, nb
467     y(i) = x(i)*s
468     ENDDO
469    
470     nbo = nb
471     ma_mult_r11 = 0
472     !-----------------------
473     END FUNCTION ma_mult_r11
474    
475     !************************************************
476    
477     INTEGER FUNCTION ma_div_r11(nb, x, s, nbo, y)
478     INTEGER :: nb, nbo
479     REAL :: x(nb), s, y(nbo)
480    
481     INTEGER :: i
482     !---------------------------------------------------------------------
483     DO i=1, nb
484     y(i) = x(i)/s
485     ENDDO
486    
487     nbo = nb
488     ma_div_r11 = 0
489     !-----------------------
490     END FUNCTION ma_div_r11
491    
492     !************************************************
493    
494     INTEGER FUNCTION ma_divi_r11(nb, x, s, nbo, y)
495     INTEGER :: nb, nbo
496     REAL :: x(nb), s, y(nbo)
497    
498     INTEGER :: i
499     !---------------------------------------------------------------------
500     DO i=1, nb
501     y(i) = s/x(i)
502     ENDDO
503    
504     nbo = nb
505     ma_divi_r11 = 0
506     !-----------------------
507     END FUNCTION ma_divi_r11
508    
509     !************************************************
510    
511     INTEGER FUNCTION ma_power_r11(nb, x, s, nbo, y)
512     INTEGER :: nb, nbo
513     REAL :: x(nb), s, y(nbo)
514    
515     INTEGER :: i
516     !---------------------------------------------------------------------
517     DO i=1, nb
518     y(i) = x(i)**s
519     ENDDO
520    
521     nbo = nb
522     ma_power_r11 = 0
523     !-----------------------
524     END FUNCTION ma_power_r11
525    
526     !************************************************
527    
528     INTEGER FUNCTION ma_fumin_r11(nb, x, s, nbo, y)
529     INTEGER :: nb, nbo
530     REAL :: x(nb), s, y(nbo)
531    
532     INTEGER :: i
533     !---------------------------------------------------------------------
534     DO i=1, nb
535     y(i) = MIN(x(i), s)
536     ENDDO
537    
538     nbo = nb
539     ma_fumin_r11 = 0
540     !------------------------
541     END FUNCTION ma_fumin_r11
542    
543     !************************************************
544    
545     INTEGER FUNCTION ma_fumax_r11(nb, x, s, nbo, y)
546     INTEGER :: nb, nbo
547     REAL :: x(nb), s, y(nbo)
548    
549     INTEGER :: i
550     !---------------------------------------------------------------------
551     DO i=1, nb
552     y(i) = MAX(x(i), s)
553     ENDDO
554    
555     nbo = nb
556     ma_fumax_r11 = 0
557     !------------------------
558     END FUNCTION ma_fumax_r11
559    
560     !************************************************
561    
562     INTEGER FUNCTION ma_fuscat_r11(nb, x, nbi, ind, miss_val, nbo, y)
563     INTEGER :: nb, nbo, nbi
564     INTEGER :: ind(nbi)
565     REAL :: x(nb), miss_val, y(nbo)
566    
567     INTEGER :: i, ii, ipos
568     !---------------------------------------------------------------------
569     ma_fuscat_r11 = 0
570    
571     y(1:nbo) = miss_val
572    
573     IF (nbi <= nb) THEN
574     ipos = 0
575     DO i=1, nbi
576     IF (ind(i) <= nbo .AND. ind(i) > 0) THEN
577     ipos = ipos+1
578     y(ind(i)) = x(ipos)
579     ELSE
580     IF (ind(i) > nbo) ma_fuscat_r11 = ma_fuscat_r11+1
581     ENDIF
582     ENDDO
583     !-- Repeat the data if needed
584     IF (MINVAL(ind) < 0) THEN
585     DO i=1, nbi
586     IF (ind(i) <= 0) THEN
587     DO ii=1, ABS(ind(i))-1
588     IF (ind(i+1)+ii <= nbo) THEN
589     y(ind(i+1)+ii) = y(ind(i+1))
590     ELSE
591     ma_fuscat_r11 = ma_fuscat_r11+1
592     ENDIF
593     ENDDO
594     ENDIF
595     ENDDO
596     ENDIF
597     ELSE
598     ma_fuscat_r11 = 1
599     ENDIF
600     !-------------------------
601     END FUNCTION ma_fuscat_r11
602    
603     !************************************************
604    
605     INTEGER FUNCTION ma_fugath_r11(nb, x, nbi, ind, miss_val, nbo, y)
606     INTEGER :: nb, nbo, nbi
607     INTEGER :: ind(nbi)
608     REAL :: x(nb), miss_val, y(nbo)
609    
610     INTEGER :: i, ipos
611     !---------------------------------------------------------------------
612     IF (nbi <= nbo) THEN
613     ma_fugath_r11 = 0
614     y(1:nbo) = miss_val
615     ipos = 0
616     DO i=1, nbi
617     IF (ipos+1 <= nbo) THEN
618     IF (ind(i) > 0) THEN
619     ipos = ipos+1
620     y(ipos) = x(ind(i))
621     ENDIF
622     ELSE
623     IF (ipos+1 > nbo) ma_fugath_r11 = ma_fugath_r11+1
624     ENDIF
625     ENDDO
626     ELSE
627     ma_fugath_r11 = 1
628     ENDIF
629    
630     nbo = ipos
631     !-------------------------
632     END FUNCTION ma_fugath_r11
633    
634     !************************************************
635    
636     INTEGER FUNCTION ma_fufill_r11(nb, x, nbi, ind, miss_val, nbo, y)
637     INTEGER :: nb, nbo, nbi
638     INTEGER :: ind(nbi)
639     REAL :: x(nb), miss_val, y(nbo)
640    
641     INTEGER :: i, ii, ipos
642     !---------------------------------------------------------------------
643     ma_fufill_r11 = 0
644    
645     IF (nbi <= nb) THEN
646     ipos = 0
647     DO i=1, nbi
648     IF (ind(i) <= nbo .AND. ind(i) > 0) THEN
649     ipos = ipos+1
650     y(ind(i)) = x(ipos)
651     ELSE
652     IF (ind(i) > nbo) ma_fufill_r11 = ma_fufill_r11+1
653     ENDIF
654     ENDDO
655     !-- Repeat the data if needed
656     IF (MINVAL(ind) < 0) THEN
657     DO i=1, nbi
658     IF (ind(i) <= 0) THEN
659     DO ii=1, ABS(ind(i))-1
660     IF (ind(i+1)+ii <= nbo) THEN
661     y(ind(i+1)+ii) = y(ind(i+1))
662     ELSE
663     ma_fufill_r11 = ma_fufill_r11+1
664     ENDIF
665     ENDDO
666     ENDIF
667     ENDDO
668     ENDIF
669     ELSE
670     ma_fufill_r11 = 1
671     ENDIF
672     !-------------------------
673     END FUNCTION ma_fufill_r11
674    
675     !************************************************
676    
677     INTEGER FUNCTION ma_fucoll_r11(nb, x, nbi, ind, miss_val, nbo, y)
678     INTEGER :: nb, nbo, nbi
679     INTEGER :: ind(nbi)
680     REAL :: x(nb), miss_val, y(nbo)
681    
682     INTEGER :: i, ipos
683     !---------------------------------------------------------------------
684     IF (nbi <= nbo) THEN
685     ma_fucoll_r11 = 0
686     ipos = 0
687     DO i=1, nbi
688     IF (ipos+1 <= nbo) THEN
689     IF (ind(i) > 0) THEN
690     ipos = ipos+1
691     y(ipos) = x(ind(i))
692     ENDIF
693     ELSE
694     IF (ipos+1 > nbo) ma_fucoll_r11 = ma_fucoll_r11+1
695     ENDIF
696     ENDDO
697     ELSE
698     ma_fucoll_r11 = 1
699     ENDIF
700    
701     nbo = ipos
702     !-------------------------
703     END FUNCTION ma_fucoll_r11
704    
705     !************************************************
706    
707     INTEGER FUNCTION ma_fuundef_r11(nb, x, nbi, ind, miss_val, nbo, y)
708     INTEGER :: nb, nbo, nbi
709     INTEGER :: ind(nbi)
710     REAL :: x(nb), miss_val, y(nbo)
711    
712     INTEGER :: i
713     !---------------------------------------------------------------------
714     IF (nbi <= nbo .AND. nbo == nb) THEN
715     ma_fuundef_r11 = 0
716     DO i=1, nbo
717     y(i) = x(i)
718     ENDDO
719     DO i=1, nbi
720     IF (ind(i) <= nbo .AND. ind(i) > 0) THEN
721     y(ind(i)) = miss_val
722     ELSE
723     IF (ind(i) > nbo) ma_fuundef_r11 = ma_fuundef_r11+1
724     ENDIF
725     ENDDO
726     ELSE
727     ma_fuundef_r11 = 1
728     ENDIF
729     !--------------------------
730     END FUNCTION ma_fuundef_r11
731    
732     !************************************************
733    
734     INTEGER FUNCTION ma_fuonly_r11(nb, x, nbi, ind, miss_val, nbo, y)
735     INTEGER :: nb, nbo, nbi
736     INTEGER :: ind(nbi)
737     REAL :: x(nb), miss_val, y(nbo)
738    
739     INTEGER :: i
740     !---------------------------------------------------------------------
741     IF ( (nbi <= nbo).AND.(nbo == nb) &
742     & .AND.ALL(ind(1:nbi) <= nbo) ) THEN
743     ma_fuonly_r11 = 0
744     y(1:nbo) = miss_val
745     DO i=1, nbi
746     IF (ind(i) > 0) THEN
747     y(ind(i)) = x(ind(i))
748     ENDIF
749     ENDDO
750     ELSE
751     ma_fuonly_r11 = 1
752     ENDIF
753     !-------------------------
754     END FUNCTION ma_fuonly_r11
755    
756     !************************************************
757    
758     !************************************************
759    
760     SUBROUTINE mathop_r21 &
761     & (fun, nb, work_in, miss_val, nb_index, nindex, scal, nb_max, work_out)
762     !- This subroutines gives an interface to the various operations
763     !- which are allowed. The interface is general enough to allow its use
764     !- for other cases.
765    
766     !- INPUT
767    
768     !- fun : function to be applied to the vector of data
769     !- nb : Length of input vector
770     !- work_in : Input vector of data (REAL)
771     !- miss_val : The value of the missing data flag (it has to be a
772     !- maximum value, in f90 : huge( a real ))
773     !- nb_index : Length of index vector
774     !- nindex : Vector of indices
775     !- scal : A scalar value for vector/scalar operations
776     !- nb_max : maximum length of output vector
777    
778     !- OUTPUT
779    
780     !- nb_max : Actual length of output variable
781     !- work_out : Output vector after the operation was applied
782     USE errioipsl, ONLY : histerr
783    
784     CHARACTER(LEN=7) :: fun
785     INTEGER :: nb(2), nb_max, nb_index
786     INTEGER :: nindex(nb_index)
787     REAL :: work_in(nb(1), nb(2)), scal, miss_val
788     REAL :: work_out(nb_max)
789    
790     INTEGER :: ierr
791    
792     INTRINSIC SIN, COS, TAN, ASIN, ACOS, ATAN, EXP, ALOG, SQRT, ABS
793     !---------------------------------------------------------------------
794     ierr = 0
795    
796     IF (scal >= miss_val-1.) THEN
797     IF (INDEX(indexfu, fun(1:LEN_TRIM(fun))) == 0) THEN
798     SELECT CASE (fun)
799     CASE('sin')
800     ierr = ma_sin_r21(nb, work_in, nb_max, work_out)
801     CASE('cos')
802     ierr = ma_cos_r21(nb, work_in, nb_max, work_out)
803     CASE('tan')
804     ierr = ma_tan_r21(nb, work_in, nb_max, work_out)
805     CASE('asin')
806     ierr = ma_asin_r21(nb, work_in, nb_max, work_out)
807     CASE('acos')
808     ierr = ma_acos_r21(nb, work_in, nb_max, work_out)
809     CASE('atan')
810     ierr = ma_atan_r21(nb, work_in, nb_max, work_out)
811     CASE('exp')
812     ierr = ma_exp_r21(nb, work_in, nb_max, work_out)
813     CASE('log')
814     ierr = ma_alog_r21(nb, work_in, nb_max, work_out)
815     CASE('sqrt')
816     ierr = ma_sqrt_r21(nb, work_in, nb_max, work_out)
817     CASE('chs')
818     ierr = ma_chs_r21(nb, work_in, nb_max, work_out)
819     CASE('abs')
820     ierr = ma_abs_r21(nb, work_in, nb_max, work_out)
821     CASE('cels')
822     ierr = ma_cels_r21(nb, work_in, nb_max, work_out)
823     CASE('kelv')
824     ierr = ma_kelv_r21(nb, work_in, nb_max, work_out)
825     CASE('deg')
826     ierr = ma_deg_r21(nb, work_in, nb_max, work_out)
827     CASE('rad')
828     ierr = ma_rad_r21(nb, work_in, nb_max, work_out)
829     CASE('ident')
830     ierr = ma_ident_r21(nb, work_in, nb_max, work_out)
831     CASE DEFAULT
832     CALL histerr(3, "mathop", &
833     & 'scalar variable undefined and no indexing', &
834     & 'but still unknown function', fun)
835     END SELECT
836     IF (ierr > 0) THEN
837     CALL histerr(3, "mathop", &
838     & 'Error while executing a simple function', fun, ' ')
839     ENDIF
840     ELSE
841     SELECT CASE (fun)
842     CASE('gather')
843     ierr = ma_fugath_r21(nb, work_in, nb_index, nindex, &
844     & miss_val, nb_max, work_out)
845     CASE('scatter')
846     IF (nb_index > (nb(1)*nb(2)) ) THEN
847     work_out(1:nb_max) = miss_val
848     ierr=1
849     ELSE
850     ierr = ma_fuscat_r21(nb, work_in, nb_index, nindex, &
851     & miss_val, nb_max, work_out)
852     ENDIF
853     CASE('coll')
854     ierr = ma_fucoll_r21(nb, work_in, nb_index, nindex, &
855     & miss_val, nb_max, work_out)
856     CASE('fill')
857     ierr = ma_fufill_r21(nb, work_in, nb_index, nindex, &
858     & miss_val, nb_max, work_out)
859     CASE('undef')
860     ierr = ma_fuundef_r21(nb, work_in, nb_index, nindex, &
861     & miss_val, nb_max, work_out)
862     CASE('only')
863     ierr = ma_fuonly_r21(nb, work_in, nb_index, nindex, &
864     & miss_val, nb_max, work_out)
865     CASE DEFAULT
866     CALL histerr(3, "mathop", &
867     & 'scalar variable undefined and indexing', &
868     & 'was requested but with unknown function', fun)
869     END SELECT
870     IF (ierr > 0) THEN
871     CALL histerr(3, "mathop_r21", &
872     & 'Error while executing an indexing function', fun, ' ')
873     ENDIF
874     ENDIF
875     ELSE
876     SELECT CASE (fun)
877     CASE('fumin')
878     ierr = ma_fumin_r21(nb, work_in, scal, nb_max, work_out)
879     CASE('fumax')
880     ierr = ma_fumax_r21(nb, work_in, scal, nb_max, work_out)
881     CASE('add')
882     ierr = ma_add_r21(nb, work_in, scal, nb_max, work_out)
883     CASE('subi')
884     ierr = ma_subi_r21(nb, work_in, scal, nb_max, work_out)
885     CASE('sub')
886     ierr = ma_sub_r21(nb, work_in, scal, nb_max, work_out)
887     CASE('mult')
888     ierr = ma_mult_r21(nb, work_in, scal, nb_max, work_out)
889     CASE('div')
890     ierr = ma_div_r21(nb, work_in, scal, nb_max, work_out)
891     CASE('divi')
892     ierr = ma_divi_r21(nb, work_in, scal, nb_max, work_out)
893     CASE('power')
894     ierr = ma_power_r21(nb, work_in, scal, nb_max, work_out)
895     CASE DEFAULT
896     CALL histerr(3, "mathop", &
897     & 'Unknown operation with a scalar', fun, ' ')
898     END SELECT
899     IF (ierr > 0) THEN
900     CALL histerr(3, "mathop", &
901     & 'Error while executing a scalar function', fun, ' ')
902     ENDIF
903     ENDIF
904     !------------------------
905     END SUBROUTINE mathop_r21
906     !-
907     !=== FUNCTIONS (only one argument)
908     !-
909     INTEGER FUNCTION ma_sin_r21(nb, x, nbo, y)
910     INTEGER :: nb(2), nbo, i, j, ij
911     REAL :: x(nb(1), nb(2)), y(nbo)
912     !---------------------------------------------------------------------
913     ij = 0
914     DO j=1, nb(2)
915     DO i=1, nb(1)
916     ij = ij+1
917     y(ij) = SIN(x(i, j))
918     ENDDO
919     ENDDO
920    
921     nbo = nb(1)*nb(2)
922     ma_sin_r21 = 0
923     !----------------------
924     END FUNCTION ma_sin_r21
925    
926     !************************************************
927    
928     INTEGER FUNCTION ma_cos_r21(nb, x, nbo, y)
929     INTEGER :: nb(2), nbo, i, j, ij
930     REAL :: x(nb(1), nb(2)), y(nbo)
931     !---------------------------------------------------------------------
932     ij = 0
933     DO j=1, nb(2)
934     DO i=1, nb(1)
935     ij = ij+1
936     y(ij) = COS(x(i, j))
937     ENDDO
938     ENDDO
939    
940     nbo = nb(1)*nb(2)
941     ma_cos_r21 = 0
942     !----------------------
943     END FUNCTION ma_cos_r21
944    
945     !************************************************
946    
947     INTEGER FUNCTION ma_tan_r21(nb, x, nbo, y)
948     INTEGER :: nb(2), nbo, i, j, ij
949     REAL :: x(nb(1), nb(2)), y(nbo)
950     !---------------------------------------------------------------------
951     ij = 0
952     DO j=1, nb(2)
953     DO i=1, nb(1)
954     ij = ij+1
955     y(ij) = TAN(x(i, j))
956     ENDDO
957     ENDDO
958    
959     nbo = nb(1)*nb(2)
960     ma_tan_r21 = 0
961     !----------------------
962     END FUNCTION ma_tan_r21
963    
964     !************************************************
965    
966     INTEGER FUNCTION ma_asin_r21(nb, x, nbo, y)
967     INTEGER :: nb(2), nbo, i, j, ij
968     REAL :: x(nb(1), nb(2)), y(nbo)
969     !---------------------------------------------------------------------
970     ij = 0
971     DO j=1, nb(2)
972     DO i=1, nb(1)
973     ij = ij+1
974     y(ij) = ASIN(x(i, j))
975     ENDDO
976     ENDDO
977    
978     nbo = nb(1)*nb(2)
979     ma_asin_r21 = 0
980     !-----------------------
981     END FUNCTION ma_asin_r21
982    
983     !************************************************
984    
985     INTEGER FUNCTION ma_acos_r21(nb, x, nbo, y)
986     INTEGER :: nb(2), nbo, i, j, ij
987     REAL :: x(nb(1), nb(2)), y(nbo)
988     !---------------------------------------------------------------------
989     ij = 0
990     DO j=1, nb(2)
991     DO i=1, nb(1)
992     ij = ij+1
993     y(ij) = ACOS(x(i, j))
994     ENDDO
995     ENDDO
996    
997     nbo = nb(1)*nb(2)
998     ma_acos_r21 = 0
999     !-----------------------
1000     END FUNCTION ma_acos_r21
1001    
1002     !************************************************
1003    
1004     INTEGER FUNCTION ma_atan_r21(nb, x, nbo, y)
1005     INTEGER :: nb(2), nbo, i, j, ij
1006     REAL :: x(nb(1), nb(2)), y(nbo)
1007     !---------------------------------------------------------------------
1008     ij = 0
1009     DO j=1, nb(2)
1010     DO i=1, nb(1)
1011     ij = ij+1
1012     y(ij) = ATAN(x(i, j))
1013     ENDDO
1014     ENDDO
1015    
1016     nbo = nb(1)*nb(2)
1017     ma_atan_r21 = 0
1018     !-----------------------
1019     END FUNCTION ma_atan_r21
1020    
1021     !************************************************
1022    
1023     INTEGER FUNCTION ma_exp_r21(nb, x, nbo, y)
1024     INTEGER :: nb(2), nbo, i, j, ij
1025     REAL :: x(nb(1), nb(2)), y(nbo)
1026     !---------------------------------------------------------------------
1027     ij = 0
1028     DO j=1, nb(2)
1029     DO i=1, nb(1)
1030     ij = ij+1
1031     y(ij) = EXP(x(i, j))
1032     ENDDO
1033     ENDDO
1034    
1035     nbo = nb(1)*nb(2)
1036     ma_exp_r21 = 0
1037     !----------------------
1038     END FUNCTION ma_exp_r21
1039    
1040     !************************************************
1041    
1042     INTEGER FUNCTION ma_alog_r21(nb, x, nbo, y)
1043     INTEGER :: nb(2), nbo, i, j, ij
1044     REAL :: x(nb(1), nb(2)), y(nbo)
1045     !---------------------------------------------------------------------
1046     ij = 0
1047     DO j=1, nb(2)
1048     DO i=1, nb(1)
1049     ij = ij+1
1050     y(ij) = ALOG(x(i, j))
1051     ENDDO
1052     ENDDO
1053    
1054     nbo = nb(1)*nb(2)
1055     ma_alog_r21 = 0
1056     !-----------------------
1057     END FUNCTION ma_alog_r21
1058    
1059     !************************************************
1060    
1061     INTEGER FUNCTION ma_sqrt_r21(nb, x, nbo, y)
1062     INTEGER :: nb(2), nbo, i, j, ij
1063     REAL :: x(nb(1), nb(2)), y(nbo)
1064     !---------------------------------------------------------------------
1065     ij = 0
1066     DO j=1, nb(2)
1067     DO i=1, nb(1)
1068     ij = ij+1
1069     y(ij) = SQRT(x(i, j))
1070     ENDDO
1071     ENDDO
1072    
1073     nbo = nb(1)*nb(2)
1074     ma_sqrt_r21 = 0
1075     !-----------------------
1076     END FUNCTION ma_sqrt_r21
1077    
1078     !************************************************
1079    
1080     INTEGER FUNCTION ma_abs_r21(nb, x, nbo, y)
1081     INTEGER :: nb(2), nbo, i, j, ij
1082     REAL :: x(nb(1), nb(2)), y(nbo)
1083     !---------------------------------------------------------------------
1084     ij = 0
1085     DO j=1, nb(2)
1086     DO i=1, nb(1)
1087     ij = ij+1
1088     y(ij) = ABS(x(i, j))
1089     ENDDO
1090     ENDDO
1091    
1092     nbo = nb(1)*nb(2)
1093     ma_abs_r21 = 0
1094     !----------------------
1095     END FUNCTION ma_abs_r21
1096    
1097     !************************************************
1098    
1099     INTEGER FUNCTION ma_chs_r21(nb, x, nbo, y)
1100     INTEGER :: nb(2), nbo, i, j, ij
1101     REAL :: x(nb(1), nb(2)), y(nbo)
1102     !---------------------------------------------------------------------
1103     ij = 0
1104     DO j=1, nb(2)
1105     DO i=1, nb(1)
1106     ij = ij+1
1107     y(ij) = x(i, j)*(-1.)
1108     ENDDO
1109     ENDDO
1110    
1111     nbo = nb(1)*nb(2)
1112     ma_chs_r21 = 0
1113     !----------------------
1114     END FUNCTION ma_chs_r21
1115    
1116     !************************************************
1117    
1118     INTEGER FUNCTION ma_cels_r21(nb, x, nbo, y)
1119     INTEGER :: nb(2), nbo, i, j, ij
1120     REAL :: x(nb(1), nb(2)), y(nbo)
1121     !---------------------------------------------------------------------
1122     ij = 0
1123     DO j=1, nb(2)
1124     DO i=1, nb(1)
1125     ij = ij+1
1126     y(ij) = x(i, j)-273.15
1127     ENDDO
1128     ENDDO
1129    
1130     nbo = nb(1)*nb(2)
1131     ma_cels_r21 = 0
1132     !-----------------------
1133     END FUNCTION ma_cels_r21
1134    
1135     !************************************************
1136    
1137     INTEGER FUNCTION ma_kelv_r21(nb, x, nbo, y)
1138     INTEGER :: nb(2), nbo, i, j, ij
1139     REAL :: x(nb(1), nb(2)), y(nbo)
1140     !---------------------------------------------------------------------
1141     ij = 0
1142     DO j=1, nb(2)
1143     DO i=1, nb(1)
1144     ij = ij+1
1145     y(ij) = x(i, j)+273.15
1146     ENDDO
1147     ENDDO
1148    
1149     nbo = nb(1)*nb(2)
1150     ma_kelv_r21 = 0
1151     !-----------------------
1152     END FUNCTION ma_kelv_r21
1153    
1154     !************************************************
1155    
1156     INTEGER FUNCTION ma_deg_r21(nb, x, nbo, y)
1157     INTEGER :: nb(2), nbo, i, j, ij
1158     REAL :: x(nb(1), nb(2)), y(nbo)
1159     !---------------------------------------------------------------------
1160     ij = 0
1161     DO j=1, nb(2)
1162     DO i=1, nb(1)
1163     ij = ij+1
1164     y(ij) = x(i, j)*57.29577951
1165     ENDDO
1166     ENDDO
1167    
1168     nbo = nb(1)*nb(2)
1169     ma_deg_r21 = 0
1170     !----------------------
1171     END FUNCTION ma_deg_r21
1172    
1173     !************************************************
1174    
1175     INTEGER FUNCTION ma_rad_r21(nb, x, nbo, y)
1176     INTEGER :: nb(2), nbo, i, j, ij
1177     REAL :: x(nb(1), nb(2)), y(nbo)
1178     !---------------------------------------------------------------------
1179     ij = 0
1180     DO j=1, nb(2)
1181     DO i=1, nb(1)
1182     ij = ij+1
1183     y(ij) = x(i, j)*0.01745329252
1184     ENDDO
1185     ENDDO
1186    
1187     nbo = nb(1)*nb(2)
1188     ma_rad_r21 = 0
1189     !----------------------
1190     END FUNCTION ma_rad_r21
1191    
1192     !************************************************
1193    
1194     INTEGER FUNCTION ma_ident_r21(nb, x, nbo, y)
1195     INTEGER :: nb(2), nbo, i, j, ij
1196     REAL :: x(nb(1), nb(2)), y(nbo)
1197     !---------------------------------------------------------------------
1198     ij = 0
1199     DO j=1, nb(2)
1200     DO i=1, nb(1)
1201     ij = ij+1
1202     y(ij) = x(i, j)
1203     ENDDO
1204     ENDDO
1205    
1206     nbo = nb(1)*nb(2)
1207     ma_ident_r21 = 0
1208     !------------------------
1209     END FUNCTION ma_ident_r21
1210     !-
1211     !=== OPERATIONS (two argument)
1212     !-
1213     INTEGER FUNCTION ma_add_r21(nb, x, s, nbo, y)
1214     INTEGER :: nb(2), nbo
1215     REAL :: x(nb(1), nb(2)), s, y(nbo)
1216    
1217     INTEGER :: i, j, ij
1218     !---------------------------------------------------------------------
1219     ij = 0
1220     DO j=1, nb(2)
1221     DO i=1, nb(1)
1222     ij = ij+1
1223     y(ij) = x(i, j)+s
1224     ENDDO
1225     ENDDO
1226    
1227     nbo = nb(1)*nb(2)
1228     ma_add_r21 = 0
1229     !----------------------
1230     END FUNCTION ma_add_r21
1231    
1232     !************************************************
1233    
1234     INTEGER FUNCTION ma_sub_r21(nb, x, s, nbo, y)
1235     INTEGER :: nb(2), nbo
1236     REAL :: x(nb(1), nb(2)), s, y(nbo)
1237    
1238     INTEGER :: i, j, ij
1239     !---------------------------------------------------------------------
1240     ij = 0
1241     DO j=1, nb(2)
1242     DO i=1, nb(1)
1243     ij = ij+1
1244     y(ij) = x(i, j)-s
1245     ENDDO
1246     ENDDO
1247    
1248     nbo = nb(1)*nb(2)
1249     ma_sub_r21 = 0
1250     !----------------------
1251     END FUNCTION ma_sub_r21
1252    
1253     !************************************************
1254    
1255     INTEGER FUNCTION ma_subi_r21(nb, x, s, nbo, y)
1256     INTEGER :: nb(2), nbo
1257     REAL :: x(nb(1), nb(2)), s, y(nbo)
1258    
1259     INTEGER :: i, j, ij
1260     !---------------------------------------------------------------------
1261     ij = 0
1262     DO j=1, nb(2)
1263     DO i=1, nb(1)
1264     ij = ij+1
1265     y(ij) = s-x(i, j)
1266     ENDDO
1267     ENDDO
1268    
1269     nbo = nb(1)*nb(2)
1270     ma_subi_r21 = 0
1271     !-----------------------
1272     END FUNCTION ma_subi_r21
1273    
1274     !************************************************
1275    
1276     INTEGER FUNCTION ma_mult_r21(nb, x, s, nbo, y)
1277     INTEGER :: nb(2), nbo
1278     REAL :: x(nb(1), nb(2)), s, y(nbo)
1279    
1280     INTEGER :: i, j, ij
1281     !---------------------------------------------------------------------
1282     ij = 0
1283     DO j=1, nb(2)
1284     DO i=1, nb(1)
1285     ij = ij+1
1286     y(ij) = x(i, j)*s
1287     ENDDO
1288     ENDDO
1289    
1290     nbo = nb(1)*nb(2)
1291     ma_mult_r21 = 0
1292     !-----------------------
1293     END FUNCTION ma_mult_r21
1294    
1295     !************************************************
1296    
1297     INTEGER FUNCTION ma_div_r21(nb, x, s, nbo, y)
1298     INTEGER :: nb(2), nbo
1299     REAL :: x(nb(1), nb(2)), s, y(nbo)
1300    
1301     INTEGER :: i, j, ij
1302     !---------------------------------------------------------------------
1303     ij = 0
1304     DO j=1, nb(2)
1305     DO i=1, nb(1)
1306     ij = ij+1
1307     y(ij) = x(i, j)/s
1308     ENDDO
1309     ENDDO
1310    
1311     nbo = nb(1)*nb(2)
1312     ma_div_r21 = 0
1313     !----------------------
1314     END FUNCTION ma_div_r21
1315    
1316     !************************************************
1317    
1318     INTEGER FUNCTION ma_divi_r21(nb, x, s, nbo, y)
1319     INTEGER :: nb(2), nbo
1320     REAL :: x(nb(1), nb(2)), s, y(nbo)
1321    
1322     INTEGER :: i, j, ij
1323     !---------------------------------------------------------------------
1324     ij = 0
1325     DO j=1, nb(2)
1326     DO i=1, nb(1)
1327     ij = ij+1
1328     y(ij) = s/x(i, j)
1329     ENDDO
1330     ENDDO
1331    
1332     nbo = nb(1)*nb(2)
1333     ma_divi_r21 = 0
1334     !-----------------------
1335     END FUNCTION ma_divi_r21
1336    
1337     !************************************************
1338    
1339     INTEGER FUNCTION ma_power_r21(nb, x, s, nbo, y)
1340     INTEGER :: nb(2), nbo
1341     REAL :: x(nb(1), nb(2)), s, y(nbo)
1342    
1343     INTEGER :: i, j, ij
1344     !---------------------------------------------------------------------
1345     ij = 0
1346     DO j=1, nb(2)
1347     DO i=1, nb(1)
1348     ij = ij+1
1349     y(ij) = x(i, j) ** s
1350     ENDDO
1351     ENDDO
1352    
1353     nbo = nb(1)*nb(2)
1354     ma_power_r21 = 0
1355     !------------------------
1356     END FUNCTION ma_power_r21
1357    
1358     !************************************************
1359    
1360     INTEGER FUNCTION ma_fumin_r21(nb, x, s, nbo, y)
1361     INTEGER :: nb(2), nbo
1362     REAL :: x(nb(1), nb(2)), s, y(nbo)
1363    
1364     INTEGER :: i, j, ij
1365     !---------------------------------------------------------------------
1366     ij = 0
1367     DO j=1, nb(2)
1368     DO i=1, nb(1)
1369     ij = ij+1
1370     y(ij) = MIN(x(i, j), s)
1371     ENDDO
1372     ENDDO
1373    
1374     nbo = nb(1)*nb(2)
1375     ma_fumin_r21 = 0
1376     !------------------------
1377     END FUNCTION ma_fumin_r21
1378    
1379     !************************************************
1380    
1381     INTEGER FUNCTION ma_fumax_r21(nb, x, s, nbo, y)
1382     INTEGER :: nb(2), nbo
1383     REAL :: x(nb(1), nb(2)), s, y(nbo)
1384    
1385     INTEGER :: i, j, ij
1386     !---------------------------------------------------------------------
1387     ij = 0
1388     DO j=1, nb(2)
1389     DO i=1, nb(1)
1390     ij = ij+1
1391     y(ij) = MAX(x(i, j), s)
1392     ENDDO
1393     ENDDO
1394    
1395     nbo = nb(1)*nb(2)
1396     ma_fumax_r21 = 0
1397     !------------------------
1398     END FUNCTION ma_fumax_r21
1399    
1400     !************************************************
1401    
1402     INTEGER FUNCTION ma_fuscat_r21(nb, x, nbi, ind, miss_val, nbo, y)
1403     INTEGER :: nb(2), nbo, nbi
1404     INTEGER :: ind(nbi)
1405     REAL :: x(nb(1), nb(2)), miss_val, y(nbo)
1406    
1407     INTEGER :: i, j, ij, ii, ipos
1408     !---------------------------------------------------------------------
1409     ma_fuscat_r21 = 0
1410    
1411     y(1:nbo) = miss_val
1412    
1413     IF (nbi <= nb(1)*nb(2)) THEN
1414     ipos = 0
1415     DO ij=1, nbi
1416     IF (ind(ij) <= nbo .AND. ind(ij) > 0) THEN
1417     ipos = ipos+1
1418     j = ((ipos-1)/nb(1))+1
1419     i = (ipos-(j-1)*nb(1))
1420     y(ind(ij)) = x(i, j)
1421     ELSE
1422     IF (ind(ij) > nbo) ma_fuscat_r21 = ma_fuscat_r21+1
1423     ENDIF
1424     ENDDO
1425     !-- Repeat the data if needed
1426     IF (MINVAL(ind) < 0) THEN
1427     DO i=1, nbi
1428     IF (ind(i) <= 0) THEN
1429     DO ii=1, ABS(ind(i))-1
1430     IF (ind(i+1)+ii <= nbo) THEN
1431     y(ind(i+1)+ii) = y(ind(i+1))
1432     ELSE
1433     ma_fuscat_r21 = ma_fuscat_r21+1
1434     ENDIF
1435     ENDDO
1436     ENDIF
1437     ENDDO
1438     ENDIF
1439     ELSE
1440     ma_fuscat_r21 = 1
1441     ENDIF
1442     !-------------------------
1443     END FUNCTION ma_fuscat_r21
1444    
1445     !************************************************
1446    
1447     INTEGER FUNCTION ma_fugath_r21(nb, x, nbi, ind, miss_val, nbo, y)
1448     INTEGER :: nb(2), nbo, nbi
1449     INTEGER :: ind(nbi)
1450     REAL :: x(nb(1), nb(2)), miss_val, y(nbo)
1451    
1452     INTEGER :: i, j, ij, ipos
1453     !---------------------------------------------------------------------
1454     IF (nbi <= nbo) THEN
1455     ma_fugath_r21 = 0
1456     y(1:nbo) = miss_val
1457     ipos = 0
1458     DO ij=1, nbi
1459     IF (ipos+1 <= nbo) THEN
1460     IF (ind(ij) > 0) THEN
1461     j = ((ind(ij)-1)/nb(1))+1
1462     i = (ind(ij)-(j-1)*nb(1))
1463     ipos = ipos+1
1464     y(ipos) = x(i, j)
1465     ENDIF
1466     ELSE
1467     IF (ipos+1 > nbo) ma_fugath_r21 = ma_fugath_r21+1
1468     ENDIF
1469     ENDDO
1470     ELSE
1471     ma_fugath_r21 = 1
1472     ENDIF
1473     nbo = ipos
1474     !-------------------------
1475     END FUNCTION ma_fugath_r21
1476    
1477     !************************************************
1478    
1479     INTEGER FUNCTION ma_fufill_r21(nb, x, nbi, ind, miss_val, nbo, y)
1480     INTEGER :: nb(2), nbo, nbi
1481     INTEGER :: ind(nbi)
1482     REAL :: x(nb(1), nb(2)), miss_val, y(nbo)
1483    
1484     INTEGER :: i, j, ij, ii, ipos
1485     !---------------------------------------------------------------------
1486     ma_fufill_r21 = 0
1487    
1488     IF (nbi <= nb(1)*nb(2)) THEN
1489     ipos = 0
1490     DO ij=1, nbi
1491     IF (ind(ij) <= nbo .AND. ind(ij) > 0) THEN
1492     ipos = ipos+1
1493     j = ((ipos-1)/nb(1))+1
1494     i = (ipos-(j-1)*nb(1))
1495     y(ind(ij)) = x(i, j)
1496     ELSE
1497     IF (ind(ij) > nbo) ma_fufill_r21 = ma_fufill_r21+1
1498     ENDIF
1499     ENDDO
1500     !-- Repeat the data if needed
1501     IF (MINVAL(ind) < 0) THEN
1502     DO i=1, nbi
1503     IF (ind(i) <= 0) THEN
1504     DO ii=1, ABS(ind(i))-1
1505     IF (ind(i+1)+ii <= nbo) THEN
1506     y(ind(i+1)+ii) = y(ind(i+1))
1507     ELSE
1508     ma_fufill_r21 = ma_fufill_r21+1
1509     ENDIF
1510     ENDDO
1511     ENDIF
1512     ENDDO
1513     ENDIF
1514     ELSE
1515     ma_fufill_r21 = 1
1516     ENDIF
1517     !-------------------------
1518     END FUNCTION ma_fufill_r21
1519    
1520     !************************************************
1521    
1522     INTEGER FUNCTION ma_fucoll_r21(nb, x, nbi, ind, miss_val, nbo, y)
1523     INTEGER :: nb(2), nbo, nbi
1524     INTEGER :: ind(nbi)
1525     REAL :: x(nb(1), nb(2)), miss_val, y(nbo)
1526    
1527     INTEGER :: i, j, ij, ipos
1528     !---------------------------------------------------------------------
1529     IF (nbi <= nbo) THEN
1530     ma_fucoll_r21 = 0
1531     ipos = 0
1532     DO ij=1, nbi
1533     IF (ipos+1 <= nbo) THEN
1534     IF (ind(ij) > 0) THEN
1535     j = ((ind(ij)-1)/nb(1))+1
1536     i = (ind(ij)-(j-1)*nb(1))
1537     ipos = ipos+1
1538     y(ipos) = x(i, j)
1539     ENDIF
1540     ELSE
1541     IF (ipos+1 > nbo) ma_fucoll_r21 = ma_fucoll_r21+1
1542     ENDIF
1543     ENDDO
1544     ELSE
1545     ma_fucoll_r21 = 1
1546     ENDIF
1547     nbo = ipos
1548     !-------------------------
1549     END FUNCTION ma_fucoll_r21
1550    
1551     !************************************************
1552    
1553     INTEGER FUNCTION ma_fuundef_r21(nb, x, nbi, ind, miss_val, nbo, y)
1554     INTEGER :: nb(2), nbo, nbi
1555     INTEGER :: ind(nbi)
1556     REAL :: x(nb(1), nb(2)), miss_val, y(nbo)
1557    
1558     INTEGER :: i, j, ij
1559     !---------------------------------------------------------------------
1560     IF (nbi <= nbo .AND. nbo == nb(1)*nb(2)) THEN
1561     ma_fuundef_r21 = 0
1562     DO ij=1, nbo
1563     j = ((ij-1)/nb(1))+1
1564     i = (ij-(j-1)*nb(1))
1565     y(ij) = x(i, j)
1566     ENDDO
1567     DO i=1, nbi
1568     IF (ind(i) <= nbo .AND. ind(i) > 0) THEN
1569     y(ind(i)) = miss_val
1570     ELSE
1571     IF (ind(i) > nbo) ma_fuundef_r21 = ma_fuundef_r21+1
1572     ENDIF
1573     ENDDO
1574     ELSE
1575     ma_fuundef_r21 = 1
1576     ENDIF
1577     !--------------------------
1578     END FUNCTION ma_fuundef_r21
1579    
1580     !************************************************
1581    
1582     INTEGER FUNCTION ma_fuonly_r21(nb, x, nbi, ind, miss_val, nbo, y)
1583     INTEGER :: nb(2), nbo, nbi
1584     INTEGER :: ind(nbi)
1585     REAL :: x(nb(1), nb(2)), miss_val, y(nbo)
1586    
1587     INTEGER :: i, j, ij
1588     !---------------------------------------------------------------------
1589     IF ( (nbi <= nbo).AND.(nbo == nb(1)*nb(2)) &
1590     & .AND.ALL(ind(1:nbi) <= nbo) ) THEN
1591     ma_fuonly_r21 = 0
1592     y(1:nbo) = miss_val
1593     DO ij=1, nbi
1594     IF (ind(ij) > 0) THEN
1595     j = ((ind(ij)-1)/nb(1))+1
1596     i = (ind(ij)-(j-1)*nb(1))
1597     y(ind(ij)) = x(i, j)
1598     ENDIF
1599     ENDDO
1600     ELSE
1601     ma_fuonly_r21 = 1
1602     ENDIF
1603     !-------------------------
1604     END FUNCTION ma_fuonly_r21
1605    
1606     !************************************************
1607    
1608     !************************************************
1609    
1610     SUBROUTINE mathop_r31 &
1611     & (fun, nb, work_in, miss_val, nb_index, nindex, scal, nb_max, work_out)
1612     !- This subroutines gives an interface to the various operations
1613     !- which are allowed. The interface is general enough to allow its use
1614     !- for other cases.
1615    
1616     !- INPUT
1617    
1618     !- fun : function to be applied to the vector of data
1619     !- nb : Length of input vector
1620     !- work_in : Input vector of data (REAL)
1621     !- miss_val : The value of the missing data flag (it has to be a
1622     !- maximum value, in f90 : huge( a real ))
1623     !- nb_index : Length of index vector
1624     !- nindex : Vector of indices
1625     !- scal : A scalar value for vector/scalar operations
1626     !- nb_max : maximum length of output vector
1627    
1628     !- OUTPUT
1629    
1630     !- nb_max : Actual length of output variable
1631     !- work_out : Output vector after the operation was applied
1632     USE errioipsl, ONLY : histerr
1633    
1634     CHARACTER(LEN=7) :: fun
1635     INTEGER :: nb(3), nb_max, nb_index
1636     INTEGER :: nindex(nb_index)
1637     REAL :: work_in(nb(1), nb(2), nb(3)), scal, miss_val
1638     REAL :: work_out(nb_max)
1639    
1640     INTEGER :: ierr
1641    
1642     INTRINSIC SIN, COS, TAN, ASIN, ACOS, ATAN, EXP, ALOG, SQRT, ABS
1643     !---------------------------------------------------------------------
1644     ierr = 0
1645    
1646     IF (scal >= miss_val-1.) THEN
1647     IF (INDEX(indexfu, fun(1:LEN_TRIM(fun))) == 0) THEN
1648     SELECT CASE (fun)
1649     CASE('sin')
1650     ierr = ma_sin_r31(nb, work_in, nb_max, work_out)
1651     CASE('cos')
1652     ierr = ma_cos_r31(nb, work_in, nb_max, work_out)
1653     CASE('tan')
1654     ierr = ma_tan_r31(nb, work_in, nb_max, work_out)
1655     CASE('asin')
1656     ierr = ma_asin_r31(nb, work_in, nb_max, work_out)
1657     CASE('acos')
1658     ierr = ma_acos_r31(nb, work_in, nb_max, work_out)
1659     CASE('atan')
1660     ierr = ma_atan_r31(nb, work_in, nb_max, work_out)
1661     CASE('exp')
1662     ierr = ma_exp_r31(nb, work_in, nb_max, work_out)
1663     CASE('log')
1664     ierr = ma_alog_r31(nb, work_in, nb_max, work_out)
1665     CASE('sqrt')
1666     ierr = ma_sqrt_r31(nb, work_in, nb_max, work_out)
1667     CASE('chs')
1668     ierr = ma_chs_r31(nb, work_in, nb_max, work_out)
1669     CASE('abs')
1670     ierr = ma_abs_r31(nb, work_in, nb_max, work_out)
1671     CASE('cels')
1672     ierr = ma_cels_r31(nb, work_in, nb_max, work_out)
1673     CASE('kelv')
1674     ierr = ma_kelv_r31(nb, work_in, nb_max, work_out)
1675     CASE('deg')
1676     ierr = ma_deg_r31(nb, work_in, nb_max, work_out)
1677     CASE('rad')
1678     ierr = ma_rad_r31(nb, work_in, nb_max, work_out)
1679     CASE('ident')
1680     ierr = ma_ident_r31(nb, work_in, nb_max, work_out)
1681     CASE DEFAULT
1682     CALL histerr(3, "mathop", &
1683     & 'scalar variable undefined and no indexing', &
1684     & 'but still unknown function', fun)
1685     END SELECT
1686     IF (ierr > 0) THEN
1687     CALL histerr(3, "mathop", &
1688     & 'Error while executing a simple function', fun, ' ')
1689     ENDIF
1690     ELSE
1691     SELECT CASE (fun)
1692     CASE('gather')
1693     ierr = ma_fugath_r31(nb, work_in, nb_index, nindex, &
1694     & miss_val, nb_max, work_out)
1695     CASE('scatter')
1696     IF (nb_index > (nb(1)*nb(2)*nb(3))) THEN
1697     work_out(1:nb_max) = miss_val
1698     ierr=1
1699     ELSE
1700     ierr = ma_fuscat_r31(nb, work_in, nb_index, nindex, &
1701     & miss_val, nb_max, work_out)
1702     ENDIF
1703     CASE('coll')
1704     ierr = ma_fucoll_r31(nb, work_in, nb_index, nindex, &
1705     & miss_val, nb_max, work_out)
1706     CASE('fill')
1707     ierr = ma_fufill_r31(nb, work_in, nb_index, nindex, &
1708     & miss_val, nb_max, work_out)
1709     CASE('undef')
1710     ierr = ma_fuundef_r31(nb, work_in, nb_index, nindex, &
1711     & miss_val, nb_max, work_out)
1712     CASE('only')
1713     ierr = ma_fuonly_r31(nb, work_in, nb_index, nindex, &
1714     & miss_val, nb_max, work_out)
1715     CASE DEFAULT
1716     CALL histerr(3, "mathop", &
1717     & 'scalar variable undefined and indexing', &
1718     & 'was requested but with unknown function', fun)
1719     END SELECT
1720     IF (ierr > 0) THEN
1721     CALL histerr(3, "mathop_r31", &
1722     & 'Error while executing an indexing function', fun, ' ')
1723     ENDIF
1724     ENDIF
1725     ELSE
1726     SELECT CASE (fun)
1727     CASE('fumin')
1728     ierr = ma_fumin_r31(nb, work_in, scal, nb_max, work_out)
1729     CASE('fumax')
1730     ierr = ma_fumax_r31(nb, work_in, scal, nb_max, work_out)
1731     CASE('add')
1732     ierr = ma_add_r31(nb, work_in, scal, nb_max, work_out)
1733     CASE('subi')
1734     ierr = ma_subi_r31(nb, work_in, scal, nb_max, work_out)
1735     CASE('sub')
1736     ierr = ma_sub_r31(nb, work_in, scal, nb_max, work_out)
1737     CASE('mult')
1738     ierr = ma_mult_r31(nb, work_in, scal, nb_max, work_out)
1739     CASE('div')
1740     ierr = ma_div_r31(nb, work_in, scal, nb_max, work_out)
1741     CASE('divi')
1742     ierr = ma_divi_r31(nb, work_in, scal, nb_max, work_out)
1743     CASE('power')
1744     ierr = ma_power_r31(nb, work_in, scal, nb_max, work_out)
1745     CASE DEFAULT
1746     CALL histerr(3, "mathop", &
1747     & 'Unknown operation with a scalar', fun, ' ')
1748     END SELECT
1749     IF (ierr > 0) THEN
1750     CALL histerr(3, "mathop", &
1751     & 'Error while executing a scalar function', fun, ' ')
1752     ENDIF
1753     ENDIF
1754     !------------------------
1755     END SUBROUTINE mathop_r31
1756     !-
1757     !=== FUNCTIONS (only one argument)
1758     !-
1759     INTEGER FUNCTION ma_sin_r31(nb, x, nbo, y)
1760     INTEGER :: nb(3), nbo, i, j, k, ij
1761     REAL :: x(nb(1), nb(2), nb(3)), y(nbo)
1762     !---------------------------------------------------------------------
1763     ij = 0
1764     DO k=1, nb(3)
1765     DO j=1, nb(2)
1766     DO i=1, nb(1)
1767     ij = ij+1
1768     y(ij) = SIN(x(i, j, k))
1769     ENDDO
1770     ENDDO
1771     ENDDO
1772    
1773     nbo = nb(1)*nb(2)*nb(3)
1774     ma_sin_r31 = 0
1775     !----------------------
1776     END FUNCTION ma_sin_r31
1777    
1778     !************************************************
1779    
1780     INTEGER FUNCTION ma_cos_r31(nb, x, nbo, y)
1781     INTEGER :: nb(3), nbo, i, j, k, ij
1782     REAL :: x(nb(1), nb(2), nb(3)), y(nbo)
1783     !---------------------------------------------------------------------
1784     ij = 0
1785     DO k=1, nb(3)
1786     DO j=1, nb(2)
1787     DO i=1, nb(1)
1788     ij = ij+1
1789     y(ij) = COS(x(i, j, k))
1790     ENDDO
1791     ENDDO
1792     ENDDO
1793    
1794     nbo = nb(1)*nb(2)*nb(3)
1795     ma_cos_r31 = 0
1796     !----------------------
1797     END FUNCTION ma_cos_r31
1798    
1799     !************************************************
1800    
1801     INTEGER FUNCTION ma_tan_r31(nb, x, nbo, y)
1802     INTEGER :: nb(3), nbo, i, j, k, ij
1803     REAL :: x(nb(1), nb(2), nb(3)), y(nbo)
1804     !---------------------------------------------------------------------
1805     ij = 0
1806     DO k=1, nb(3)
1807     DO j=1, nb(2)
1808     DO i=1, nb(1)
1809     ij = ij+1
1810     y(ij) = TAN(x(i, j, k))
1811     ENDDO
1812     ENDDO
1813     ENDDO
1814    
1815     nbo = nb(1)*nb(2)*nb(3)
1816     ma_tan_r31 = 0
1817     !----------------------
1818     END FUNCTION ma_tan_r31
1819    
1820     !************************************************
1821    
1822     INTEGER FUNCTION ma_asin_r31(nb, x, nbo, y)
1823     INTEGER :: nb(3), nbo, i, j, k, ij
1824     REAL :: x(nb(1), nb(2), nb(3)), y(nbo)
1825     !---------------------------------------------------------------------
1826     ij = 0
1827     DO k=1, nb(3)
1828     DO j=1, nb(2)
1829     DO i=1, nb(1)
1830     ij = ij+1
1831     y(ij) = ASIN(x(i, j, k))
1832     ENDDO
1833     ENDDO
1834     ENDDO
1835    
1836     nbo = nb(1)*nb(2)*nb(3)
1837     ma_asin_r31 = 0
1838     !-----------------------
1839     END FUNCTION ma_asin_r31
1840    
1841     !************************************************
1842    
1843     INTEGER FUNCTION ma_acos_r31(nb, x, nbo, y)
1844     INTEGER :: nb(3), nbo, i, j, k, ij
1845     REAL :: x(nb(1), nb(2), nb(3)), y(nbo)
1846     !---------------------------------------------------------------------
1847     ij = 0
1848     DO k=1, nb(3)
1849     DO j=1, nb(2)
1850     DO i=1, nb(1)
1851     ij = ij+1
1852     y(ij) = ACOS(x(i, j, k))
1853     ENDDO
1854     ENDDO
1855     ENDDO
1856    
1857     nbo = nb(1)*nb(2)*nb(3)
1858     ma_acos_r31 = 0
1859     !-----------------------
1860     END FUNCTION ma_acos_r31
1861    
1862     !************************************************
1863    
1864     INTEGER FUNCTION ma_atan_r31(nb, x, nbo, y)
1865     INTEGER :: nb(3), nbo, i, j, k, ij
1866     REAL :: x(nb(1), nb(2), nb(3)), y(nbo)
1867     !---------------------------------------------------------------------
1868     ij = 0
1869     DO k=1, nb(3)
1870     DO j=1, nb(2)
1871     DO i=1, nb(1)
1872     ij = ij+1
1873     y(ij) = ATAN(x(i, j, k))
1874     ENDDO
1875     ENDDO
1876     ENDDO
1877    
1878     nbo = nb(1)*nb(2)*nb(3)
1879     ma_atan_r31 = 0
1880     !-----------------------
1881     END FUNCTION ma_atan_r31
1882    
1883     !************************************************
1884    
1885     INTEGER FUNCTION ma_exp_r31(nb, x, nbo, y)
1886     INTEGER :: nb(3), nbo, i, j, k, ij
1887     REAL :: x(nb(1), nb(2), nb(3)), y(nbo)
1888     !---------------------------------------------------------------------
1889     ij = 0
1890     DO k=1, nb(3)
1891     DO j=1, nb(2)
1892     DO i=1, nb(1)
1893     ij = ij+1
1894     y(ij) = EXP(x(i, j, k))
1895     ENDDO
1896     ENDDO
1897     ENDDO
1898    
1899     nbo = nb(1)*nb(2)*nb(3)
1900     ma_exp_r31 = 0
1901     !----------------------
1902     END FUNCTION ma_exp_r31
1903    
1904     !************************************************
1905    
1906     INTEGER FUNCTION ma_alog_r31(nb, x, nbo, y)
1907     INTEGER :: nb(3), nbo, i, j, k, ij
1908     REAL :: x(nb(1), nb(2), nb(3)), y(nbo)
1909     !---------------------------------------------------------------------
1910     ij = 0
1911     DO k=1, nb(3)
1912     DO j=1, nb(2)
1913     DO i=1, nb(1)
1914     ij = ij+1
1915     y(ij) = ALOG(x(i, j, k))
1916     ENDDO
1917     ENDDO
1918     ENDDO
1919    
1920     nbo = nb(1)*nb(2)*nb(3)
1921     ma_alog_r31 = 0
1922     !-----------------------
1923     END FUNCTION ma_alog_r31
1924    
1925     !************************************************
1926    
1927     INTEGER FUNCTION ma_sqrt_r31(nb, x, nbo, y)
1928     INTEGER :: nb(3), nbo, i, j, k, ij
1929     REAL :: x(nb(1), nb(2), nb(3)), y(nbo)
1930     !---------------------------------------------------------------------
1931     ij = 0
1932     DO k=1, nb(3)
1933     DO j=1, nb(2)
1934     DO i=1, nb(1)
1935     ij = ij+1
1936     y(ij) = SQRT(x(i, j, k))
1937     ENDDO
1938     ENDDO
1939     ENDDO
1940    
1941     nbo = nb(1)*nb(2)*nb(3)
1942     ma_sqrt_r31 = 0
1943     !-----------------------
1944     END FUNCTION ma_sqrt_r31
1945    
1946     !************************************************
1947    
1948     INTEGER FUNCTION ma_abs_r31(nb, x, nbo, y)
1949     INTEGER :: nb(3), nbo, i, j, k, ij
1950     REAL :: x(nb(1), nb(2), nb(3)), y(nbo)
1951     !---------------------------------------------------------------------
1952     ij = 0
1953     DO k=1, nb(3)
1954     DO j=1, nb(2)
1955     DO i=1, nb(1)
1956     ij = ij+1
1957     y(ij) = ABS(x(i, j, k))
1958     ENDDO
1959     ENDDO
1960     ENDDO
1961    
1962     nbo = nb(1)*nb(2)*nb(3)
1963     ma_abs_r31 = 0
1964     !----------------------
1965     END FUNCTION ma_abs_r31
1966    
1967     !************************************************
1968    
1969     INTEGER FUNCTION ma_chs_r31(nb, x, nbo, y)
1970     INTEGER :: nb(3), nbo, i, j, k, ij
1971     REAL :: x(nb(1), nb(2), nb(3)), y(nbo)
1972     !---------------------------------------------------------------------
1973     ij = 0
1974     DO k=1, nb(3)
1975     DO j=1, nb(2)
1976     DO i=1, nb(1)
1977     ij = ij+1
1978     y(ij) = x(i, j, k)*(-1.)
1979     ENDDO
1980     ENDDO
1981     ENDDO
1982    
1983     nbo = nb(1)*nb(2)*nb(3)
1984     ma_chs_r31 = 0
1985     !----------------------
1986     END FUNCTION ma_chs_r31
1987    
1988     !************************************************
1989    
1990     INTEGER FUNCTION ma_cels_r31(nb, x, nbo, y)
1991     INTEGER :: nb(3), nbo, i, j, k, ij
1992     REAL :: x(nb(1), nb(2), nb(3)), y(nbo)
1993     !---------------------------------------------------------------------
1994     ij = 0
1995     DO k=1, nb(3)
1996     DO j=1, nb(2)
1997     DO i=1, nb(1)
1998     ij = ij+1
1999     y(ij) = x(i, j, k)-273.15
2000     ENDDO
2001     ENDDO
2002     ENDDO
2003    
2004     nbo = nb(1)*nb(2)*nb(3)
2005     ma_cels_r31 = 0
2006     !-----------------------
2007     END FUNCTION ma_cels_r31
2008    
2009     !************************************************
2010    
2011     INTEGER FUNCTION ma_kelv_r31(nb, x, nbo, y)
2012     INTEGER :: nb(3), nbo, i, j, k, ij
2013     REAL :: x(nb(1), nb(2), nb(3)), y(nbo)
2014     !---------------------------------------------------------------------
2015     ij = 0
2016     DO k=1, nb(3)
2017     DO j=1, nb(2)
2018     DO i=1, nb(1)
2019     ij = ij+1
2020     y(ij) = x(i, j, k)+273.15
2021     ENDDO
2022     ENDDO
2023     ENDDO
2024    
2025     nbo = nb(1)*nb(2)*nb(3)
2026     ma_kelv_r31 = 0
2027     !-----------------------
2028     END FUNCTION ma_kelv_r31
2029    
2030     !************************************************
2031    
2032     INTEGER FUNCTION ma_deg_r31(nb, x, nbo, y)
2033     INTEGER :: nb(3), nbo, i, j, k, ij
2034     REAL :: x(nb(1), nb(2), nb(3)), y(nbo)
2035     !---------------------------------------------------------------------
2036     ij = 0
2037     DO k=1, nb(3)
2038     DO j=1, nb(2)
2039     DO i=1, nb(1)
2040     ij = ij+1
2041     y(ij) = x(i, j, k)*57.29577951
2042     ENDDO
2043     ENDDO
2044     ENDDO
2045    
2046     nbo = nb(1)*nb(2)*nb(3)
2047     ma_deg_r31 = 0
2048     !----------------------
2049     END FUNCTION ma_deg_r31
2050    
2051     !************************************************
2052    
2053     INTEGER FUNCTION ma_rad_r31(nb, x, nbo, y)
2054     INTEGER :: nb(3), nbo, i, j, k, ij
2055     REAL :: x(nb(1), nb(2), nb(3)), y(nbo)
2056     !---------------------------------------------------------------------
2057     ij = 0
2058     DO k=1, nb(3)
2059     DO j=1, nb(2)
2060     DO i=1, nb(1)
2061     ij = ij+1
2062     y(ij) = x(i, j, k)*0.01745329252
2063     ENDDO
2064     ENDDO
2065     ENDDO
2066    
2067     nbo = nb(1)*nb(2)*nb(3)
2068     ma_rad_r31 = 0
2069     !----------------------
2070     END FUNCTION ma_rad_r31
2071    
2072     !************************************************
2073    
2074     INTEGER FUNCTION ma_ident_r31(nb, x, nbo, y)
2075     INTEGER :: nb(3), nbo, i, j, k, ij
2076     REAL :: x(nb(1), nb(2), nb(3)), y(nbo)
2077     !---------------------------------------------------------------------
2078     ij = 0
2079     DO k=1, nb(3)
2080     DO j=1, nb(2)
2081     DO i=1, nb(1)
2082     ij = ij+1
2083     y(ij) = x(i, j, k)
2084     ENDDO
2085     ENDDO
2086     ENDDO
2087    
2088     nbo = nb(1)*nb(2)*nb(3)
2089     ma_ident_r31 = 0
2090     !------------------------
2091     END FUNCTION ma_ident_r31
2092     !-
2093     !=== OPERATIONS (two argument)
2094     !-
2095     INTEGER FUNCTION ma_add_r31(nb, x, s, nbo, y)
2096     INTEGER :: nb(3), nbo
2097     REAL :: x(nb(1), nb(2), nb(3)), s, y(nbo)
2098    
2099     INTEGER :: i, j, k, ij
2100     !---------------------------------------------------------------------
2101     ij = 0
2102     DO k=1, nb(3)
2103     DO j=1, nb(2)
2104     DO i=1, nb(1)
2105     ij = ij+1
2106     y(ij) = x(i, j, k)+s
2107     ENDDO
2108     ENDDO
2109     ENDDO
2110    
2111     nbo = nb(1)*nb(2)*nb(3)
2112     ma_add_r31 = 0
2113     !----------------------
2114     END FUNCTION ma_add_r31
2115    
2116     !************************************************
2117    
2118     INTEGER FUNCTION ma_sub_r31(nb, x, s, nbo, y)
2119     INTEGER :: nb(3), nbo
2120     REAL :: x(nb(1), nb(2), nb(3)), s, y(nbo)
2121    
2122     INTEGER :: i, j, k, ij
2123     !---------------------------------------------------------------------
2124     ij = 0
2125     DO k=1, nb(3)
2126     DO j=1, nb(2)
2127     DO i=1, nb(1)
2128     ij = ij+1
2129     y(ij) = x(i, j, k)-s
2130     ENDDO
2131     ENDDO
2132     ENDDO
2133    
2134     nbo = nb(1)*nb(2)*nb(3)
2135     ma_sub_r31 = 0
2136     !----------------------
2137     END FUNCTION ma_sub_r31
2138    
2139     !************************************************
2140    
2141     INTEGER FUNCTION ma_subi_r31(nb, x, s, nbo, y)
2142     INTEGER :: nb(3), nbo
2143     REAL :: x(nb(1), nb(2), nb(3)), s, y(nbo)
2144    
2145     INTEGER :: i, j, k, ij
2146     !---------------------------------------------------------------------
2147     ij = 0
2148     DO k=1, nb(3)
2149     DO j=1, nb(2)
2150     DO i=1, nb(1)
2151     ij = ij+1
2152     y(ij) = s-x(i, j, k)
2153     ENDDO
2154     ENDDO
2155     ENDDO
2156    
2157     nbo = nb(1)*nb(2)*nb(3)
2158     ma_subi_r31 = 0
2159     !-----------------------
2160     END FUNCTION ma_subi_r31
2161    
2162     !************************************************
2163    
2164     INTEGER FUNCTION ma_mult_r31(nb, x, s, nbo, y)
2165     INTEGER :: nb(3), nbo
2166     REAL :: x(nb(1), nb(2), nb(3)), s, y(nbo)
2167    
2168     INTEGER :: i, j, k, ij
2169     !---------------------------------------------------------------------
2170     ij = 0
2171     DO k=1, nb(3)
2172     DO j=1, nb(2)
2173     DO i=1, nb(1)
2174     ij = ij+1
2175     y(ij) = x(i, j, k)*s
2176     ENDDO
2177     ENDDO
2178     ENDDO
2179    
2180     nbo = nb(1)*nb(2)*nb(3)
2181     ma_mult_r31 = 0
2182     !-----------------------
2183     END FUNCTION ma_mult_r31
2184    
2185     !************************************************
2186    
2187     INTEGER FUNCTION ma_div_r31(nb, x, s, nbo, y)
2188     INTEGER :: nb(3), nbo
2189     REAL :: x(nb(1), nb(2), nb(3)), s, y(nbo)
2190    
2191     INTEGER :: i, j, k, ij
2192     !---------------------------------------------------------------------
2193     ij = 0
2194     DO k=1, nb(3)
2195     DO j=1, nb(2)
2196     DO i=1, nb(1)
2197     ij = ij+1
2198     y(ij) = x(i, j, k)/s
2199     ENDDO
2200     ENDDO
2201     ENDDO
2202    
2203     nbo = nb(1)*nb(2)*nb(3)
2204     ma_div_r31 = 0
2205     !----------------------
2206     END FUNCTION ma_div_r31
2207    
2208     !************************************************
2209    
2210     INTEGER FUNCTION ma_divi_r31(nb, x, s, nbo, y)
2211     INTEGER :: nb(3), nbo
2212     REAL :: x(nb(1), nb(2), nb(3)), s, y(nbo)
2213    
2214     INTEGER :: i, j, k, ij
2215     !---------------------------------------------------------------------
2216     ij = 0
2217     DO k=1, nb(3)
2218     DO j=1, nb(2)
2219     DO i=1, nb(1)
2220     ij = ij+1
2221     y(ij) = s/x(i, j, k)
2222     ENDDO
2223     ENDDO
2224     ENDDO
2225    
2226     nbo = nb(1)*nb(2)*nb(3)
2227     ma_divi_r31 = 0
2228     !-----------------------
2229     END FUNCTION ma_divi_r31
2230    
2231     !************************************************
2232    
2233     INTEGER FUNCTION ma_power_r31(nb, x, s, nbo, y)
2234     INTEGER :: nb(3), nbo
2235     REAL :: x(nb(1), nb(2), nb(3)), s, y(nbo)
2236    
2237     INTEGER :: i, j, k, ij
2238     !---------------------------------------------------------------------
2239     ij = 0
2240     DO k=1, nb(3)
2241     DO j=1, nb(2)
2242     DO i=1, nb(1)
2243     ij = ij+1
2244     y(ij) = x(i, j, k)**s
2245     ENDDO
2246     ENDDO
2247     ENDDO
2248    
2249     nbo = nb(1)*nb(2)*nb(3)
2250     ma_power_r31 = 0
2251     !------------------------
2252     END FUNCTION ma_power_r31
2253    
2254     !************************************************
2255    
2256     INTEGER FUNCTION ma_fumin_r31(nb, x, s, nbo, y)
2257     INTEGER :: nb(3), nbo
2258     REAL :: x(nb(1), nb(2), nb(3)), s, y(nbo)
2259    
2260     INTEGER :: i, j, k, ij
2261     !---------------------------------------------------------------------
2262     ij = 0
2263     DO k=1, nb(3)
2264     DO j=1, nb(2)
2265     DO i=1, nb(1)
2266     ij = ij+1
2267     y(ij) = MIN(x(i, j, k), s)
2268     ENDDO
2269     ENDDO
2270     ENDDO
2271    
2272     nbo = nb(1)*nb(2)*nb(3)
2273     ma_fumin_r31 = 0
2274     !------------------------
2275     END FUNCTION ma_fumin_r31
2276    
2277     !************************************************
2278    
2279     INTEGER FUNCTION ma_fumax_r31(nb, x, s, nbo, y)
2280     INTEGER :: nb(3), nbo
2281     REAL :: x(nb(1), nb(2), nb(3)), s, y(nbo)
2282    
2283     INTEGER :: i, j, k, ij
2284     !---------------------------------------------------------------------
2285     ij = 0
2286     DO k=1, nb(3)
2287     DO j=1, nb(2)
2288     DO i=1, nb(1)
2289     ij = ij+1
2290     y(ij) = MAX(x(i, j, k), s)
2291     ENDDO
2292     ENDDO
2293     ENDDO
2294    
2295     nbo = nb(1)*nb(2)*nb(3)
2296     ma_fumax_r31 = 0
2297     !------------------------
2298     END FUNCTION ma_fumax_r31
2299    
2300     !************************************************
2301    
2302     INTEGER FUNCTION ma_fuscat_r31(nb, x, nbi, ind, miss_val, nbo, y)
2303     INTEGER :: nb(3), nbo, nbi
2304     INTEGER :: ind(nbi)
2305     REAL :: x(nb(1), nb(2), nb(3)), miss_val, y(nbo)
2306    
2307     INTEGER :: i, j, k, ij, ii, ipos, ipp, isb
2308     !---------------------------------------------------------------------
2309     ma_fuscat_r31 = 0
2310    
2311     y(1:nbo) = miss_val
2312    
2313     IF (nbi <= nb(1)*nb(2)*nb(3)) THEN
2314     ipos = 0
2315     isb = nb(1)*nb(2)
2316     DO ij=1, nbi
2317     IF (ind(ij) <= nbo .AND. ind(ij) > 0) THEN
2318     ipos = ipos+1
2319     k = ((ipos-1)/isb)+1
2320     ipp = ipos-(k-1)*isb
2321     j = ((ipp-1)/nb(1))+1
2322     i = (ipp-(j-1)*nb(1))
2323     y(ind(ij)) = x(i, j, k)
2324     ELSE
2325     IF (ind(ij) > nbo) ma_fuscat_r31 = ma_fuscat_r31+1
2326     ENDIF
2327     ENDDO
2328     !-- Repeat the data if needed
2329     IF (MINVAL(ind) < 0) THEN
2330     DO i=1, nbi
2331     IF (ind(i) <= 0) THEN
2332     DO ii=1, ABS(ind(i))-1
2333     IF (ind(i+1)+ii <= nbo) THEN
2334     y(ind(i+1)+ii) = y(ind(i+1))
2335     ELSE
2336     ma_fuscat_r31 = ma_fuscat_r31+1
2337     ENDIF
2338     ENDDO
2339     ENDIF
2340     ENDDO
2341     ENDIF
2342     ELSE
2343     ma_fuscat_r31 = 1
2344     ENDIF
2345     !-------------------------
2346     END FUNCTION ma_fuscat_r31
2347    
2348     !************************************************
2349    
2350     INTEGER FUNCTION ma_fugath_r31(nb, x, nbi, ind, miss_val, nbo, y)
2351     INTEGER :: nb(3), nbo, nbi
2352     INTEGER :: ind(nbi)
2353     REAL :: x(nb(1), nb(2), nb(3)), miss_val, y(nbo)
2354    
2355     INTEGER :: i, j, k, ij, ipos, ipp, isb
2356     !---------------------------------------------------------------------
2357     IF (nbi <= nbo) THEN
2358     ma_fugath_r31 = 0
2359     y(1:nbo) = miss_val
2360     ipos = 0
2361     isb = nb(1)*nb(2)
2362     DO ij=1, nbi
2363     IF (ipos+1 <= nbo) THEN
2364     IF (ind(ij) > 0) THEN
2365     k = ((ind(ij)-1)/isb)+1
2366     ipp = ind(ij)-(k-1)*isb
2367     j = ((ipp-1)/nb(1))+1
2368     i = (ipp-(j-1)*nb(1))
2369     ipos = ipos+1
2370     y(ipos) = x(i, j, k)
2371     ENDIF
2372     ELSE
2373     IF (ipos+1 > nbo) ma_fugath_r31 = ma_fugath_r31+1
2374     ENDIF
2375     ENDDO
2376     ELSE
2377     ma_fugath_r31 = 1
2378     ENDIF
2379     nbo = ipos
2380     !-------------------------
2381     END FUNCTION ma_fugath_r31
2382    
2383     !************************************************
2384    
2385     INTEGER FUNCTION ma_fufill_r31(nb, x, nbi, ind, miss_val, nbo, y)
2386     INTEGER :: nb(3), nbo, nbi
2387     INTEGER :: ind(nbi)
2388     REAL :: x(nb(1), nb(2), nb(3)), miss_val, y(nbo)
2389    
2390     INTEGER :: i, j, k, ij, ii, ipos, ipp, isb
2391     !---------------------------------------------------------------------
2392     ma_fufill_r31 = 0
2393     IF (nbi <= nb(1)*nb(2)*nb(3)) THEN
2394     ipos = 0
2395     isb = nb(1)*nb(2)
2396     DO ij=1, nbi
2397     IF (ind(ij) <= nbo .AND. ind(ij) > 0) THEN
2398     ipos = ipos+1
2399     k = ((ipos-1)/isb)+1
2400     ipp = ipos-(k-1)*isb
2401     j = ((ipp-1)/nb(1))+1
2402     i = (ipp-(j-1)*nb(1))
2403     y(ind(ij)) = x(i, j, k)
2404     ELSE
2405     IF (ind(ij) > nbo) ma_fufill_r31 = ma_fufill_r31+1
2406     ENDIF
2407     ENDDO
2408     !-- Repeat the data if needed
2409     IF (MINVAL(ind) < 0) THEN
2410     DO i=1, nbi
2411     IF (ind(i) <= 0) THEN
2412     DO ii=1, ABS(ind(i))-1
2413     IF (ind(i+1)+ii <= nbo) THEN
2414     y(ind(i+1)+ii) = y(ind(i+1))
2415     ELSE
2416     ma_fufill_r31 = ma_fufill_r31+1
2417     ENDIF
2418     ENDDO
2419     ENDIF
2420     ENDDO
2421     ENDIF
2422     ELSE
2423     ma_fufill_r31 = 1
2424     ENDIF
2425     !-------------------------
2426     END FUNCTION ma_fufill_r31
2427    
2428     !************************************************
2429    
2430     INTEGER FUNCTION ma_fucoll_r31(nb, x, nbi, ind, miss_val, nbo, y)
2431     INTEGER :: nb(3), nbo, nbi
2432     INTEGER :: ind(nbi)
2433     REAL :: x(nb(1), nb(2), nb(3)), miss_val, y(nbo)
2434    
2435     INTEGER :: i, j, k, ij, ipos, ipp, isb
2436     !---------------------------------------------------------------------
2437     IF (nbi <= nbo) THEN
2438     ma_fucoll_r31 = 0
2439     ipos = 0
2440     isb = nb(1)*nb(2)
2441     DO ij=1, nbi
2442     IF (ipos+1 <= nbo) THEN
2443     IF (ind(ij) > 0) THEN
2444     k = ((ind(ij)-1)/isb)+1
2445     ipp = ind(ij)-(k-1)*isb
2446     j = ((ipp-1)/nb(1))+1
2447     i = (ipp-(j-1)*nb(1))
2448     ipos = ipos+1
2449     y(ipos) = x(i, j, k)
2450     ENDIF
2451     ELSE
2452     IF (ipos+1 > nbo) ma_fucoll_r31 = ma_fucoll_r31+1
2453     ENDIF
2454     ENDDO
2455     ELSE
2456     ma_fucoll_r31 = 1
2457     ENDIF
2458     nbo = ipos
2459     !-------------------------
2460     END FUNCTION ma_fucoll_r31
2461    
2462     !************************************************
2463    
2464     INTEGER FUNCTION ma_fuundef_r31(nb, x, nbi, ind, miss_val, nbo, y)
2465     INTEGER :: nb(3), nbo, nbi
2466     INTEGER :: ind(nbi)
2467     REAL :: x(nb(1), nb(2), nb(3)), miss_val, y(nbo)
2468    
2469     INTEGER :: i, j, k, ij, ipp, isb
2470     !---------------------------------------------------------------------
2471     IF (nbi <= nbo .AND. nbo == nb(1)*nb(2)*nb(3)) THEN
2472     ma_fuundef_r31 = 0
2473     isb = nb(1)*nb(2)
2474     DO ij=1, nbo
2475     k = ((ij-1)/isb)+1
2476     ipp = ij-(k-1)*isb
2477     j = ((ipp-1)/nb(1))+1
2478     i = (ipp-(j-1)*nb(1))
2479     y(ij) = x(i, j, k)
2480     ENDDO
2481     DO i=1, nbi
2482     IF (ind(i) <= nbo .AND. ind(i) > 0) THEN
2483     y(ind(i)) = miss_val
2484     ELSE
2485     IF (ind(i) > nbo) ma_fuundef_r31 = ma_fuundef_r31+1
2486     ENDIF
2487     ENDDO
2488     ELSE
2489     ma_fuundef_r31 = 1
2490     ENDIF
2491     !--------------------------
2492     END FUNCTION ma_fuundef_r31
2493    
2494     !************************************************
2495    
2496     INTEGER FUNCTION ma_fuonly_r31(nb, x, nbi, ind, miss_val, nbo, y)
2497     INTEGER :: nb(3), nbo, nbi
2498     INTEGER :: ind(nbi)
2499     REAL :: x(nb(1), nb(2), nb(3)), miss_val, y(nbo)
2500    
2501     INTEGER :: i, j, k, ij, ipp, isb
2502     !---------------------------------------------------------------------
2503     IF ( (nbi <= nbo).AND.(nbo == nb(1)*nb(2)*nb(3)) &
2504     & .AND.ALL(ind(1:nbi) <= nbo) ) THEN
2505     ma_fuonly_r31 = 0
2506     y(1:nbo) = miss_val
2507     isb = nb(1)*nb(2)
2508     DO ij=1, nbi
2509     IF (ind(ij) > 0) THEN
2510     k = ((ind(ij)-1)/isb)+1
2511     ipp = ind(ij)-(k-1)*isb
2512     j = ((ipp-1)/nb(1))+1
2513     i = (ipp-(j-1)*nb(1))
2514     y(ind(ij)) = x(i, j, k)
2515     ENDIF
2516     ENDDO
2517     ELSE
2518     ma_fuonly_r31 = 1
2519     ENDIF
2520     !-------------------------
2521     END FUNCTION ma_fuonly_r31
2522    
2523     END MODULE mathop_m
  ViewVC Help
Powered by ViewVC 1.1.21