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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 61 by guez, Tue Apr 6 17:52:58 2010 UTC revision 62 by guez, Thu Jul 26 14:37:37 2012 UTC
# Line 9  MODULE mathop_m Line 9  MODULE mathop_m
9    
10    INTERFACE mathop    INTERFACE mathop
11       MODULE PROCEDURE mathop_r11, mathop_r21, mathop_r31       MODULE PROCEDURE mathop_r11, mathop_r21, mathop_r31
12    END INTERFACE    END INTERFACE mathop
13    
14    !- Variables used to detect and identify the operations    ! Variables used to detect and identify the operations
15    CHARACTER(LEN=120), SAVE :: indexfu = &    CHARACTER(LEN=120):: indexfu = 'gather, scatter, fill, coll, undef, only'
        'gather, scatter, fill, coll, undef, only'  
16    
17  CONTAINS  CONTAINS
18    
19    SUBROUTINE mathop_r11 &    SUBROUTINE mathop_r11(fun, nb, work_in, miss_val, nb_index, nindex, scal, &
20         &  (fun, nb, work_in, miss_val, nb_index, nindex, scal, nb_max, work_out)         nb_max, work_out)
     !- This subroutines gives an interface to the various operation  
     !- which are allowed. The interface is general enough to allow its use  
     !- for other cases.  
   
     !- INPUT  
   
     !- fun      : function to be applied to the vector of data  
     !- nb       : Length of input vector  
     !- work_in  : Input vector of data (REAL)  
     !- miss_val : The value of the missing data flag (it has to be a  
     !-            maximum value, in f90 : huge( a real ))  
     !- nb_index : Length of index vector  
     !- nindex   : Vector of indices  
     !- scal     : A scalar value for vector/scalar operations  
     !- nb_max   : maximum length of output vector  
21    
22      !- OUTPUT      ! This subroutine gives an interface to the various operations
23        ! which are allowed. The interface is general enough to allow its
24        ! use 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      !- nb_max   : Actual length of output variable      ! OUTPUT
39      !- work_out : Output vector after the operation was applied  
40        ! nb_max   : Actual length of output variable
41        ! work_out : Output vector after the operation was applied
42      USE errioipsl, ONLY : histerr      USE errioipsl, ONLY : histerr
43        USE mathop2, ONLY : ma_abs_r11, ma_acos_r11, ma_add_r11, ma_alog_r11, &
44             ma_asin_r11, ma_atan_r11, ma_cels_r11, ma_chs_r11, ma_cos_r11, &
45             ma_deg_r11, ma_divi_r11, ma_div_r11, ma_exp_r11, ma_fucoll_r11, &
46             ma_fufill_r11, ma_fugath_r11, ma_fumax_r11, ma_fumin_r11, &
47             ma_fuonly_r11, ma_fuscat_r11, ma_fuundef_r11, ma_ident_r11, &
48             ma_kelv_r11, ma_mult_r11, ma_power_r11, ma_rad_r11, ma_sin_r11, &
49             ma_sqrt_r11, ma_subi_r11, ma_sub_r11, ma_tan_r11
50    
51      CHARACTER(LEN=7) :: fun      CHARACTER(LEN=7) :: fun
52      INTEGER :: nb, nb_max, nb_index      INTEGER :: nb, nb_max, nb_index
# Line 50  CONTAINS Line 57  CONTAINS
57      INTEGER :: ierr      INTEGER :: ierr
58    
59      INTRINSIC SIN, COS, TAN, ASIN, ACOS, ATAN, EXP, ALOG, SQRT, ABS      INTRINSIC SIN, COS, TAN, ASIN, ACOS, ATAN, EXP, ALOG, SQRT, ABS
60      !---------------------------------------------------------------------      !--------------------------------------------------------------------
61      ierr = 0      ierr = 0
62    
63      IF (scal >= miss_val-1.) THEN      IF (scal >= miss_val-1.) THEN
# Line 90  CONTAINS Line 97  CONTAINS
97               ierr = ma_ident_r11(nb, work_in, nb_max, work_out)               ierr = ma_ident_r11(nb, work_in, nb_max, work_out)
98            CASE DEFAULT            CASE DEFAULT
99               CALL histerr(3, "mathop", &               CALL histerr(3, "mathop", &
100                    &        'scalar variable undefined and no indexing', &          'scalar variable undefined and no indexing', &
101                    &        'but still unknown function', fun)          'but still unknown function', fun)
102            END SELECT            END SELECT
103            IF (ierr > 0) THEN            IF (ierr > 0) THEN
104               CALL histerr(3, "mathop", &               CALL histerr(3, "mathop", &
105                    &        'Error while executing a simple function', fun, ' ')          'Error while executing a simple function', fun, ' ')
106            ENDIF            ENDIF
107         ELSE         ELSE
108            SELECT CASE (fun)            SELECT CASE (fun)
109            CASE('gather')            CASE('gather')
110               ierr = ma_fugath_r11(nb, work_in, nb_index, nindex, &               ierr = ma_fugath_r11(nb, work_in, nb_index, nindex, &
111                    &                            miss_val, nb_max, work_out)                              miss_val, nb_max, work_out)
112            CASE('scatter')            CASE('scatter')
113               IF (nb_index > nb) THEN               IF (nb_index > nb) THEN
114                  work_out(1:nb_max) = miss_val                  work_out(1:nb_max) = miss_val
115                  ierr=1                  ierr=1
116               ELSE               ELSE
117                  ierr = ma_fuscat_r11(nb, work_in, nb_index, nindex, &                  ierr = ma_fuscat_r11(nb, work_in, nb_index, nindex, &
118                       &                              miss_val, nb_max, work_out)                                miss_val, nb_max, work_out)
119               ENDIF               ENDIF
120            CASE('coll')            CASE('coll')
121               ierr = ma_fucoll_r11(nb, work_in, nb_index, nindex, &               ierr = ma_fucoll_r11(nb, work_in, nb_index, nindex, &
122                    &                            miss_val, nb_max, work_out)                              miss_val, nb_max, work_out)
123            CASE('fill')            CASE('fill')
124               ierr = ma_fufill_r11(nb, work_in, nb_index, nindex, &               ierr = ma_fufill_r11(nb, work_in, nb_index, nindex, &
125                    &                            miss_val, nb_max, work_out)                              miss_val, nb_max, work_out)
126            CASE('undef')            CASE('undef')
127               ierr = ma_fuundef_r11(nb, work_in, nb_index, nindex, &               ierr = ma_fuundef_r11(nb, work_in, nb_index, nindex, &
128                    &                            miss_val, nb_max, work_out)                              miss_val, nb_max, work_out)
129            CASE('only')            CASE('only')
130               ierr = ma_fuonly_r11(nb, work_in, nb_index, nindex, &               ierr = ma_fuonly_r11(nb, work_in, nb_index, nindex, &
131                    &                            miss_val, nb_max, work_out)                              miss_val, nb_max, work_out)
132            CASE DEFAULT            CASE DEFAULT
133               CALL histerr(3, "mathop", &               CALL histerr(3, "mathop", &
134                    &        'scalar variable undefined and indexing', &          'scalar variable undefined and indexing', &
135                    &        'was requested but with unknown function', fun)          'was requested but with unknown function', fun)
136            END SELECT            END SELECT
137            IF (ierr > 0) THEN            IF (ierr > 0) THEN
138               CALL histerr(3, "mathop_r11", &               CALL histerr(3, "mathop_r11", &
139                    &        'Error while executing an indexing function', fun, ' ')          'Error while executing an indexing function', fun, ' ')
140            ENDIF            ENDIF
141         ENDIF         ENDIF
142      ELSE      ELSE
# Line 154  CONTAINS Line 161  CONTAINS
161            ierr = ma_power_r11(nb, work_in, scal, nb_max, work_out)            ierr = ma_power_r11(nb, work_in, scal, nb_max, work_out)
162         CASE DEFAULT         CASE DEFAULT
163            CALL histerr(3, "mathop", &            CALL histerr(3, "mathop", &
164                 &      'Unknown operation with a scalar', fun, ' ')        'Unknown operation with a scalar', fun, ' ')
165         END SELECT         END SELECT
166         IF (ierr > 0) THEN         IF (ierr > 0) THEN
167            CALL histerr(3, "mathop", &            CALL histerr(3, "mathop", &
168                 &      'Error while executing a scalar function', fun, ' ')        'Error while executing a scalar function', fun, ' ')
169         ENDIF         ENDIF
170      ENDIF      ENDIF
     !------------------------  
   END SUBROUTINE mathop_r11  
   !-  
   !=== FUNCTIONS (only one argument)  
   !-  
   INTEGER FUNCTION ma_sin_r11(nb, x, nbo, y)  
     INTEGER :: nb, nbo, i  
     REAL :: x(nb), y(nbo)  
     !---------------------------------------------------------------------  
     DO i=1, nb  
        y(i) = SIN(x(i))  
     ENDDO  
   
     nbo = nb  
     ma_sin_r11 = 0  
     !----------------------  
   END FUNCTION ma_sin_r11  
   
   !************************************************  
   
   INTEGER FUNCTION ma_cos_r11(nb, x, nbo, y)  
     INTEGER :: nb, nbo, i  
     REAL :: x(nb), y(nbo)  
     !---------------------------------------------------------------------  
     DO i=1, nb  
        y(i) = COS(x(i))  
     ENDDO  
   
     nbo = nb  
     ma_cos_r11 = 0  
     !----------------------  
   END FUNCTION ma_cos_r11  
   
   !************************************************  
   
   INTEGER FUNCTION ma_tan_r11(nb, x, nbo, y)  
     INTEGER :: nb, nbo, i  
     REAL :: x(nb), y(nbo)  
     !---------------------------------------------------------------------  
     DO i=1, nb  
        y(i) = TAN(x(i))  
     ENDDO  
   
     nbo = nb  
     ma_tan_r11 = 0  
     !----------------------  
   END FUNCTION ma_tan_r11  
   
   !************************************************  
   
   INTEGER FUNCTION ma_asin_r11(nb, x, nbo, y)  
     INTEGER :: nb, nbo, i  
     REAL :: x(nb), y(nbo)  
     !---------------------------------------------------------------------  
     DO i=1, nb  
        y(i) = ASIN(x(i))  
     ENDDO  
   
     nbo = nb  
     ma_asin_r11 = 0  
     !-----------------------  
   END FUNCTION ma_asin_r11  
   
   !************************************************  
   
   INTEGER FUNCTION ma_acos_r11(nb, x, nbo, y)  
     INTEGER :: nb, nbo, i  
     REAL :: x(nb), y(nbo)  
     !---------------------------------------------------------------------  
     DO i=1, nb  
        y(i) = ACOS(x(i))  
     ENDDO  
   
     nbo = nb  
     ma_acos_r11 = 0  
     !-----------------------  
   END FUNCTION ma_acos_r11  
   
   !************************************************  
   
   INTEGER FUNCTION ma_atan_r11(nb, x, nbo, y)  
     INTEGER :: nb, nbo, i  
     REAL :: x(nb), y(nbo)  
     !---------------------------------------------------------------------  
     DO i=1, nb  
        y(i) = ATAN(x(i))  
     ENDDO  
   
     nbo = nb  
     ma_atan_r11 = 0  
     !-----------------------  
   END FUNCTION ma_atan_r11  
   
   !************************************************  
   
   INTEGER FUNCTION ma_exp_r11(nb, x, nbo, y)  
     INTEGER :: nb, nbo, i  
     REAL :: x(nb), y(nbo)  
     !---------------------------------------------------------------------  
     DO i=1, nb  
        y(i) = EXP(x(i))  
     ENDDO  
   
     nbo = nb  
     ma_exp_r11 = 0  
     !----------------------  
   END FUNCTION ma_exp_r11  
   
   !************************************************  
   
   INTEGER FUNCTION ma_alog_r11(nb, x, nbo, y)  
     INTEGER :: nb, nbo, i  
     REAL :: x(nb), y(nbo)  
     !---------------------------------------------------------------------  
     DO i=1, nb  
        y(i) = alog(x(i))  
     ENDDO  
   
     nbo = nb  
     ma_alog_r11 = 0  
     !-----------------------  
   END FUNCTION ma_alog_r11  
   
   !************************************************  
   
   INTEGER FUNCTION ma_sqrt_r11(nb, x, nbo, y)  
     INTEGER :: nb, nbo, i  
     REAL :: x(nb), y(nbo)  
     !---------------------------------------------------------------------  
     DO i=1, nb  
        y(i) = SQRT(x(i))  
     ENDDO  
   
     nbo = nb  
     ma_sqrt_r11 = 0  
     !-----------------------  
   END FUNCTION ma_sqrt_r11  
   
   !************************************************  
   
   INTEGER FUNCTION ma_abs_r11(nb, x, nbo, y)  
     INTEGER :: nb, nbo, i  
     REAL :: x(nb), y(nbo)  
     !---------------------------------------------------------------------  
     DO i=1, nb  
        y(i) = ABS(x(i))  
     ENDDO  
   
     nbo = nb  
     ma_abs_r11 = 0  
     !----------------------  
   END FUNCTION ma_abs_r11  
   
   !************************************************  
   
   INTEGER FUNCTION ma_chs_r11(nb, x, nbo, y)  
     INTEGER :: nb, nbo, i  
     REAL :: x(nb), y(nbo)  
     !---------------------------------------------------------------------  
     DO i=1, nb  
        y(i) = x(i)*(-1.)  
     ENDDO  
   
     nbo = nb  
     ma_chs_r11 = 0  
     !----------------------  
   END FUNCTION ma_chs_r11  
   
   !************************************************  
   
   INTEGER FUNCTION ma_cels_r11(nb, x, nbo, y)  
     INTEGER :: nb, nbo, i  
     REAL :: x(nb), y(nbo)  
     !---------------------------------------------------------------------  
     DO i=1, nb  
        y(i) = x(i)-273.15  
     ENDDO  
   
     nbo = nb  
     ma_cels_r11 = 0  
     !-----------------------  
   END FUNCTION ma_cels_r11  
   
   !************************************************  
   
   INTEGER FUNCTION ma_kelv_r11(nb, x, nbo, y)  
     INTEGER :: nb, nbo, i  
     REAL :: x(nb), y(nbo)  
     !---------------------------------------------------------------------  
     DO i=1, nb  
        y(i) = x(i)+273.15  
     ENDDO  
   
     nbo = nb  
     ma_kelv_r11 = 0  
     !-----------------------  
   END FUNCTION ma_kelv_r11  
   
   !************************************************  
   
   INTEGER FUNCTION ma_deg_r11(nb, x, nbo, y)  
     INTEGER :: nb, nbo, i  
     REAL :: x(nb), y(nbo)  
     !---------------------------------------------------------------------  
     DO i=1, nb  
        y(i) = x(i)*57.29577951  
     ENDDO  
   
     nbo = nb  
     ma_deg_r11 = 0  
     !-----------------------  
   END FUNCTION ma_deg_r11  
   
   !************************************************  
   
   INTEGER FUNCTION ma_rad_r11(nb, x, nbo, y)  
     INTEGER :: nb, nbo, i  
     REAL :: x(nb), y(nbo)  
     !---------------------------------------------------------------------  
     DO i=1, nb  
        y(i) = x(i)*0.01745329252  
     ENDDO  
   
     nbo = nb  
     ma_rad_r11 = 0  
     !----------------------  
   END FUNCTION ma_rad_r11  
   
   !************************************************  
   
   INTEGER FUNCTION ma_ident_r11(nb, x, nbo, y)  
     INTEGER :: nb, nbo, i  
     REAL :: x(nb), y(nbo)  
     !---------------------------------------------------------------------  
     DO i=1, nb  
        y(i) = x(i)  
     ENDDO  
   
     nbo = nb  
     ma_ident_r11 = 0  
     !------------------------  
   END FUNCTION ma_ident_r11  
   !-  
   !=== OPERATIONS (two argument)  
   !-  
   INTEGER FUNCTION ma_add_r11(nb, x, s, nbo, y)  
     INTEGER :: nb, nbo  
     REAL :: x(nb), s, y(nbo)  
   
     INTEGER :: i  
     !---------------------------------------------------------------------  
     DO i=1, nb  
        y(i) = x(i)+s  
     ENDDO  
   
     nbo = nb  
     ma_add_r11 = 0  
     !-----------------------  
   END FUNCTION ma_add_r11  
   
   !************************************************  
   
   INTEGER FUNCTION ma_sub_r11(nb, x, s, nbo, y)  
     INTEGER :: nb, nbo  
     REAL :: x(nb), s, y(nbo)  
   
     INTEGER :: i  
     !---------------------------------------------------------------------  
     DO i=1, nb  
        y(i) = x(i)-s  
     ENDDO  
   
     nbo = nb  
     ma_sub_r11 = 0  
     !----------------------  
   END FUNCTION ma_sub_r11  
   
   !************************************************  
   
   INTEGER FUNCTION ma_subi_r11(nb, x, s, nbo, y)  
     INTEGER :: nb, nbo  
     REAL :: x(nb), s, y(nbo)  
   
     INTEGER :: i  
     !---------------------------------------------------------------------  
     DO i=1, nb  
        y(i) =  s-x(i)  
     ENDDO  
   
     nbo = nb  
     ma_subi_r11 = 0  
     !-----------------------  
   END FUNCTION ma_subi_r11  
   
   !************************************************  
   
   INTEGER FUNCTION ma_mult_r11(nb, x, s, nbo, y)  
     INTEGER :: nb, nbo  
     REAL :: x(nb), s, y(nbo)  
   
     INTEGER :: i  
     !---------------------------------------------------------------------  
     DO i=1, nb  
        y(i) = x(i)*s  
     ENDDO  
   
     nbo = nb  
     ma_mult_r11 = 0  
     !-----------------------  
   END FUNCTION ma_mult_r11  
   
   !************************************************  
   
   INTEGER FUNCTION ma_div_r11(nb, x, s, nbo, y)  
     INTEGER :: nb, nbo  
     REAL :: x(nb), s, y(nbo)  
   
     INTEGER :: i  
     !---------------------------------------------------------------------  
     DO i=1, nb  
        y(i) = x(i)/s  
     ENDDO  
   
     nbo = nb  
     ma_div_r11 = 0  
     !-----------------------  
   END FUNCTION ma_div_r11  
   
   !************************************************  
   
   INTEGER FUNCTION ma_divi_r11(nb, x, s, nbo, y)  
     INTEGER :: nb, nbo  
     REAL :: x(nb), s, y(nbo)  
   
     INTEGER :: i  
     !---------------------------------------------------------------------  
     DO i=1, nb  
        y(i) = s/x(i)  
     ENDDO  
   
     nbo = nb  
     ma_divi_r11 = 0  
171      !-----------------------      !-----------------------
172    END FUNCTION ma_divi_r11    END SUBROUTINE mathop_r11
173      !
174    !************************************************    SUBROUTINE mathop_r21(fun, nb, work_in, miss_val, nb_index, nindex, scal, &
175           nb_max, work_out)
176    INTEGER FUNCTION ma_power_r11(nb, x, s, nbo, y)  
177      INTEGER :: nb, nbo      ! This subroutines gives an interface to the various operations
178      REAL :: x(nb), s, y(nbo)      ! which are allowed. The interface is general enough to allow its use
179        ! for other cases.
180      INTEGER :: i  
181      !---------------------------------------------------------------------      ! INPUT
182      DO i=1, nb  
183         y(i) = x(i)**s      ! fun      : function to be applied to the vector of data
184      ENDDO      ! nb       : Length of input vector
185        ! work_in  : Input vector of data (REAL)
186      nbo = nb      ! miss_val : The value of the missing data flag (it has to be a
187      ma_power_r11 = 0      !            maximum value, in f90 : huge( a real ))
188      !-----------------------      ! nb_index : Length of index vector
189    END FUNCTION ma_power_r11      ! nindex   : Vector of indices
190        ! scal     : A scalar value for vector/scalar operations
191    !************************************************      ! nb_max   : maximum length of output vector
   
   INTEGER FUNCTION ma_fumin_r11(nb, x, s, nbo, y)  
     INTEGER :: nb, nbo  
     REAL :: x(nb), s, y(nbo)  
   
     INTEGER :: i  
     !---------------------------------------------------------------------  
     DO i=1, nb  
        y(i) = MIN(x(i), s)  
     ENDDO  
   
     nbo = nb  
     ma_fumin_r11 = 0  
     !------------------------  
   END FUNCTION ma_fumin_r11  
   
   !************************************************  
   
   INTEGER FUNCTION ma_fumax_r11(nb, x, s, nbo, y)  
     INTEGER :: nb, nbo  
     REAL :: x(nb), s, y(nbo)  
   
     INTEGER :: i  
     !---------------------------------------------------------------------  
     DO i=1, nb  
        y(i) = MAX(x(i), s)  
     ENDDO  
   
     nbo = nb  
     ma_fumax_r11 = 0  
     !------------------------  
   END FUNCTION ma_fumax_r11  
   
   !************************************************  
   
   INTEGER FUNCTION ma_fuscat_r11(nb, x, nbi, ind, miss_val, nbo, y)  
     INTEGER :: nb, nbo, nbi  
     INTEGER :: ind(nbi)  
     REAL :: x(nb), miss_val, y(nbo)  
   
     INTEGER :: i, ii, ipos  
     !---------------------------------------------------------------------  
     ma_fuscat_r11 = 0  
   
     y(1:nbo) = miss_val  
   
     IF (nbi <= nb) THEN  
        ipos = 0  
        DO i=1, nbi  
           IF (ind(i) <= nbo .AND. ind(i) > 0) THEN  
              ipos = ipos+1  
              y(ind(i)) = x(ipos)  
           ELSE  
              IF (ind(i) > nbo) ma_fuscat_r11  = ma_fuscat_r11+1  
           ENDIF  
        ENDDO  
        !-- Repeat the data if needed  
        IF (MINVAL(ind) < 0) THEN  
           DO i=1, nbi  
              IF (ind(i) <= 0) THEN  
                 DO ii=1, ABS(ind(i))-1  
                    IF (ind(i+1)+ii <= nbo) THEN  
                       y(ind(i+1)+ii) = y(ind(i+1))  
                    ELSE  
                       ma_fuscat_r11  = ma_fuscat_r11+1  
                    ENDIF  
                 ENDDO  
              ENDIF  
           ENDDO  
        ENDIF  
     ELSE  
        ma_fuscat_r11  = 1  
     ENDIF  
     !-------------------------  
   END FUNCTION ma_fuscat_r11  
   
   !************************************************  
   
   INTEGER FUNCTION ma_fugath_r11(nb, x, nbi, ind, miss_val, nbo, y)  
     INTEGER :: nb, nbo, nbi  
     INTEGER :: ind(nbi)  
     REAL :: x(nb), miss_val, y(nbo)  
   
     INTEGER :: i, ipos  
     !---------------------------------------------------------------------  
     IF (nbi <= nbo) THEN  
        ma_fugath_r11 = 0  
        y(1:nbo) = miss_val  
        ipos = 0  
        DO i=1, nbi  
           IF (ipos+1 <= nbo) THEN  
              IF (ind(i) > 0) THEN  
                 ipos = ipos+1  
                 y(ipos) = x(ind(i))  
              ENDIF  
           ELSE  
              IF (ipos+1 > nbo) ma_fugath_r11  = ma_fugath_r11+1  
           ENDIF  
        ENDDO  
     ELSE  
        ma_fugath_r11 = 1  
     ENDIF  
   
     nbo = ipos  
     !-------------------------  
   END FUNCTION ma_fugath_r11  
   
   !************************************************  
   
   INTEGER FUNCTION ma_fufill_r11(nb, x, nbi, ind, miss_val, nbo, y)  
     INTEGER :: nb, nbo, nbi  
     INTEGER :: ind(nbi)  
     REAL :: x(nb), miss_val, y(nbo)  
   
     INTEGER :: i, ii, ipos  
     !---------------------------------------------------------------------  
     ma_fufill_r11 = 0  
   
     IF (nbi <= nb) THEN  
        ipos = 0  
        DO i=1, nbi  
           IF (ind(i) <= nbo .AND. ind(i) > 0) THEN  
              ipos = ipos+1  
              y(ind(i)) = x(ipos)  
           ELSE  
              IF (ind(i) > nbo) ma_fufill_r11  = ma_fufill_r11+1  
           ENDIF  
        ENDDO  
        !-- Repeat the data if needed  
        IF (MINVAL(ind) < 0) THEN  
           DO i=1, nbi  
              IF (ind(i) <= 0) THEN  
                 DO ii=1, ABS(ind(i))-1  
                    IF (ind(i+1)+ii <= nbo) THEN  
                       y(ind(i+1)+ii) = y(ind(i+1))  
                    ELSE  
                       ma_fufill_r11  = ma_fufill_r11+1  
                    ENDIF  
                 ENDDO  
              ENDIF  
           ENDDO  
        ENDIF  
     ELSE  
        ma_fufill_r11  = 1  
     ENDIF  
     !-------------------------  
   END FUNCTION ma_fufill_r11  
   
   !************************************************  
   
   INTEGER FUNCTION ma_fucoll_r11(nb, x, nbi, ind, miss_val, nbo, y)  
     INTEGER :: nb, nbo, nbi  
     INTEGER :: ind(nbi)  
     REAL :: x(nb), miss_val, y(nbo)  
   
     INTEGER :: i, ipos  
     !---------------------------------------------------------------------  
     IF (nbi <= nbo) THEN  
        ma_fucoll_r11 = 0  
        ipos = 0  
        DO i=1, nbi  
           IF (ipos+1 <= nbo) THEN  
              IF (ind(i) > 0) THEN  
                 ipos = ipos+1  
                 y(ipos) = x(ind(i))  
              ENDIF  
           ELSE  
              IF (ipos+1 > nbo) ma_fucoll_r11  = ma_fucoll_r11+1  
           ENDIF  
        ENDDO  
     ELSE  
        ma_fucoll_r11 = 1  
     ENDIF  
   
     nbo = ipos  
     !-------------------------  
   END FUNCTION ma_fucoll_r11  
   
   !************************************************  
   
   INTEGER FUNCTION ma_fuundef_r11(nb, x, nbi, ind, miss_val, nbo, y)  
     INTEGER :: nb, nbo, nbi  
     INTEGER :: ind(nbi)  
     REAL :: x(nb), miss_val, y(nbo)  
   
     INTEGER :: i  
     !---------------------------------------------------------------------  
     IF (nbi <= nbo .AND. nbo == nb) THEN  
        ma_fuundef_r11 = 0  
        DO i=1, nbo  
           y(i) = x(i)  
        ENDDO  
        DO i=1, nbi  
           IF (ind(i) <= nbo .AND. ind(i) > 0) THEN  
              y(ind(i)) =  miss_val  
           ELSE  
              IF (ind(i) > nbo) ma_fuundef_r11  = ma_fuundef_r11+1  
           ENDIF  
        ENDDO  
     ELSE  
        ma_fuundef_r11 = 1  
     ENDIF  
     !--------------------------  
   END FUNCTION ma_fuundef_r11  
   
   !************************************************  
   
   INTEGER FUNCTION ma_fuonly_r11(nb, x, nbi, ind, miss_val, nbo, y)  
     INTEGER :: nb, nbo, nbi  
     INTEGER :: ind(nbi)  
     REAL :: x(nb), miss_val, y(nbo)  
   
     INTEGER :: i  
     !---------------------------------------------------------------------  
     IF (     (nbi <= nbo).AND.(nbo == nb) &  
          &    .AND.ALL(ind(1:nbi) <= nbo) ) THEN  
        ma_fuonly_r11 = 0  
        y(1:nbo) = miss_val  
        DO i=1, nbi  
           IF (ind(i) > 0) THEN  
              y(ind(i)) =  x(ind(i))  
           ENDIF  
        ENDDO  
     ELSE  
        ma_fuonly_r11 = 1  
     ENDIF  
     !-------------------------  
   END FUNCTION ma_fuonly_r11  
   
   !************************************************  
   
   !************************************************  
   
   SUBROUTINE mathop_r21 &  
        &  (fun, nb, work_in, miss_val, nb_index, nindex, scal, nb_max, work_out)  
     !- This subroutines gives an interface to the various operations  
     !- which are allowed. The interface is general enough to allow its use  
     !- for other cases.  
   
     !- INPUT  
   
     !- fun      : function to be applied to the vector of data  
     !- nb       : Length of input vector  
     !- work_in  : Input vector of data (REAL)  
     !- miss_val : The value of the missing data flag (it has to be a  
     !-            maximum value, in f90 : huge( a real ))  
     !- nb_index : Length of index vector  
     !- nindex   : Vector of indices  
     !- scal     : A scalar value for vector/scalar operations  
     !- nb_max   : maximum length of output vector  
192    
193      !- OUTPUT      ! OUTPUT
194    
195      !- nb_max   : Actual length of output variable      ! nb_max   : Actual length of output variable
196      !- work_out : Output vector after the operation was applied      ! work_out : Output vector after the operation was applied
197      USE errioipsl, ONLY : histerr      USE errioipsl, ONLY : histerr
198        USE mathop2, ONLY : ma_abs_r21, ma_acos_r21, ma_add_r21, ma_alog_r21, &
199             ma_asin_r21, ma_atan_r21, ma_cels_r21, ma_chs_r21, ma_cos_r21, &
200             ma_deg_r21, ma_divi_r21, ma_div_r21, ma_exp_r21, ma_fucoll_r21, &
201             ma_fufill_r21, ma_fugath_r21, ma_fumax_r21, ma_fumin_r21, &
202             ma_fuonly_r21, ma_fuscat_r21, ma_fuundef_r21, ma_ident_r21, &
203             ma_kelv_r21, ma_mult_r21, ma_power_r21, ma_rad_r21, ma_sin_r21, &
204             ma_sqrt_r21, ma_subi_r21, ma_sub_r21, ma_tan_r21
205    
206      CHARACTER(LEN=7) :: fun      CHARACTER(LEN=7) :: fun
207      INTEGER :: nb(2), nb_max, nb_index      INTEGER :: nb(2), nb_max, nb_index
# Line 790  CONTAINS Line 212  CONTAINS
212      INTEGER :: ierr      INTEGER :: ierr
213    
214      INTRINSIC SIN, COS, TAN, ASIN, ACOS, ATAN, EXP, ALOG, SQRT, ABS      INTRINSIC SIN, COS, TAN, ASIN, ACOS, ATAN, EXP, ALOG, SQRT, ABS
215      !---------------------------------------------------------------------      !--------------------------------------------------------------------
216      ierr = 0      ierr = 0
217    
218      IF (scal >= miss_val-1.) THEN      IF (scal >= miss_val-1.) THEN
# Line 830  CONTAINS Line 252  CONTAINS
252               ierr = ma_ident_r21(nb, work_in, nb_max, work_out)               ierr = ma_ident_r21(nb, work_in, nb_max, work_out)
253            CASE DEFAULT            CASE DEFAULT
254               CALL histerr(3, "mathop", &               CALL histerr(3, "mathop", &
255                    &        'scalar variable undefined and no indexing', &          'scalar variable undefined and no indexing', &
256                    &        'but still unknown function', fun)          'but still unknown function', fun)
257            END SELECT            END SELECT
258            IF (ierr > 0) THEN            IF (ierr > 0) THEN
259               CALL histerr(3, "mathop", &               CALL histerr(3, "mathop", &
260                    &        'Error while executing a simple function', fun, ' ')          'Error while executing a simple function', fun, ' ')
261            ENDIF            ENDIF
262         ELSE         ELSE
263            SELECT CASE (fun)            SELECT CASE (fun)
264            CASE('gather')            CASE('gather')
265               ierr = ma_fugath_r21(nb, work_in, nb_index, nindex, &               ierr = ma_fugath_r21(nb, work_in, nb_index, nindex, &
266                    &                            miss_val, nb_max, work_out)                              miss_val, nb_max, work_out)
267            CASE('scatter')            CASE('scatter')
268               IF (nb_index > (nb(1)*nb(2)) ) THEN               IF (nb_index > (nb(1)*nb(2)) ) THEN
269                  work_out(1:nb_max) = miss_val                  work_out(1:nb_max) = miss_val
270                  ierr=1                  ierr=1
271               ELSE               ELSE
272                  ierr = ma_fuscat_r21(nb, work_in, nb_index, nindex, &                  ierr = ma_fuscat_r21(nb, work_in, nb_index, nindex, &
273                       &                             miss_val, nb_max, work_out)                               miss_val, nb_max, work_out)
274               ENDIF               ENDIF
275            CASE('coll')            CASE('coll')
276               ierr = ma_fucoll_r21(nb, work_in, nb_index, nindex, &               ierr = ma_fucoll_r21(nb, work_in, nb_index, nindex, &
277                    &                           miss_val, nb_max, work_out)                             miss_val, nb_max, work_out)
278            CASE('fill')            CASE('fill')
279               ierr = ma_fufill_r21(nb, work_in, nb_index, nindex, &               ierr = ma_fufill_r21(nb, work_in, nb_index, nindex, &
280                    &                           miss_val, nb_max, work_out)                             miss_val, nb_max, work_out)
281            CASE('undef')            CASE('undef')
282               ierr = ma_fuundef_r21(nb, work_in, nb_index, nindex, &               ierr = ma_fuundef_r21(nb, work_in, nb_index, nindex, &
283                    &                           miss_val, nb_max, work_out)                             miss_val, nb_max, work_out)
284            CASE('only')            CASE('only')
285               ierr = ma_fuonly_r21(nb, work_in, nb_index, nindex, &               ierr = ma_fuonly_r21(nb, work_in, nb_index, nindex, &
286                    &                           miss_val, nb_max, work_out)                             miss_val, nb_max, work_out)
287            CASE DEFAULT            CASE DEFAULT
288               CALL histerr(3, "mathop", &               CALL histerr(3, "mathop", &
289                    &        'scalar variable undefined and indexing', &          'scalar variable undefined and indexing', &
290                    &        'was requested but with unknown function', fun)          'was requested but with unknown function', fun)
291            END SELECT            END SELECT
292            IF (ierr > 0) THEN            IF (ierr > 0) THEN
293               CALL histerr(3, "mathop_r21", &               CALL histerr(3, "mathop_r21", &
294                    &        'Error while executing an indexing function', fun, ' ')          'Error while executing an indexing function', fun, ' ')
295            ENDIF            ENDIF
296         ENDIF         ENDIF
297      ELSE      ELSE
# Line 894  CONTAINS Line 316  CONTAINS
316            ierr = ma_power_r21(nb, work_in, scal, nb_max, work_out)            ierr = ma_power_r21(nb, work_in, scal, nb_max, work_out)
317         CASE DEFAULT         CASE DEFAULT
318            CALL histerr(3, "mathop", &            CALL histerr(3, "mathop", &
319                 &      'Unknown operation with a scalar', fun, ' ')        'Unknown operation with a scalar', fun, ' ')
320         END SELECT         END SELECT
321         IF (ierr > 0) THEN         IF (ierr > 0) THEN
322            CALL histerr(3, "mathop", &            CALL histerr(3, "mathop", &
323                 &      'Error while executing a scalar function', fun, ' ')        'Error while executing a scalar function', fun, ' ')
324         ENDIF         ENDIF
325      ENDIF      ENDIF
     !------------------------  
   END SUBROUTINE mathop_r21  
   !-  
   !=== FUNCTIONS (only one argument)  
   !-  
   INTEGER FUNCTION ma_sin_r21(nb, x, nbo, y)  
     INTEGER :: nb(2), nbo, i, j, ij  
     REAL :: x(nb(1), nb(2)), y(nbo)  
     !---------------------------------------------------------------------  
     ij = 0  
     DO j=1, nb(2)  
        DO i=1, nb(1)  
           ij = ij+1  
           y(ij) = SIN(x(i, j))  
        ENDDO  
     ENDDO  
   
     nbo = nb(1)*nb(2)  
     ma_sin_r21 = 0  
     !----------------------  
   END FUNCTION ma_sin_r21  
   
   !************************************************  
   
   INTEGER FUNCTION ma_cos_r21(nb, x, nbo, y)  
     INTEGER :: nb(2), nbo, i, j, ij  
     REAL :: x(nb(1), nb(2)), y(nbo)  
     !---------------------------------------------------------------------  
     ij = 0  
     DO j=1, nb(2)  
        DO i=1, nb(1)  
           ij = ij+1  
           y(ij) = COS(x(i, j))  
        ENDDO  
     ENDDO  
   
     nbo = nb(1)*nb(2)  
     ma_cos_r21 = 0  
     !----------------------  
   END FUNCTION ma_cos_r21  
   
   !************************************************  
   
   INTEGER FUNCTION ma_tan_r21(nb, x, nbo, y)  
     INTEGER :: nb(2), nbo, i, j, ij  
     REAL :: x(nb(1), nb(2)), y(nbo)  
     !---------------------------------------------------------------------  
     ij = 0  
     DO j=1, nb(2)  
        DO i=1, nb(1)  
           ij = ij+1  
           y(ij) = TAN(x(i, j))  
        ENDDO  
     ENDDO  
   
     nbo = nb(1)*nb(2)  
     ma_tan_r21 = 0  
     !----------------------  
   END FUNCTION ma_tan_r21  
   
   !************************************************  
   
   INTEGER FUNCTION ma_asin_r21(nb, x, nbo, y)  
     INTEGER :: nb(2), nbo, i, j, ij  
     REAL :: x(nb(1), nb(2)), y(nbo)  
     !---------------------------------------------------------------------  
     ij = 0  
     DO j=1, nb(2)  
        DO i=1, nb(1)  
           ij = ij+1  
           y(ij) = ASIN(x(i, j))  
        ENDDO  
     ENDDO  
   
     nbo = nb(1)*nb(2)  
     ma_asin_r21 = 0  
     !-----------------------  
   END FUNCTION ma_asin_r21  
   
   !************************************************  
   
   INTEGER FUNCTION ma_acos_r21(nb, x, nbo, y)  
     INTEGER :: nb(2), nbo, i, j, ij  
     REAL :: x(nb(1), nb(2)), y(nbo)  
     !---------------------------------------------------------------------  
     ij = 0  
     DO j=1, nb(2)  
        DO i=1, nb(1)  
           ij = ij+1  
           y(ij) = ACOS(x(i, j))  
        ENDDO  
     ENDDO  
   
     nbo = nb(1)*nb(2)  
     ma_acos_r21 = 0  
     !-----------------------  
   END FUNCTION ma_acos_r21  
   
   !************************************************  
   
   INTEGER FUNCTION ma_atan_r21(nb, x, nbo, y)  
     INTEGER :: nb(2), nbo, i, j, ij  
     REAL :: x(nb(1), nb(2)), y(nbo)  
     !---------------------------------------------------------------------  
     ij = 0  
     DO j=1, nb(2)  
        DO i=1, nb(1)  
           ij = ij+1  
           y(ij) = ATAN(x(i, j))  
        ENDDO  
     ENDDO  
   
     nbo = nb(1)*nb(2)  
     ma_atan_r21 = 0  
     !-----------------------  
   END FUNCTION ma_atan_r21  
   
   !************************************************  
   
   INTEGER FUNCTION ma_exp_r21(nb, x, nbo, y)  
     INTEGER :: nb(2), nbo, i, j, ij  
     REAL :: x(nb(1), nb(2)), y(nbo)  
     !---------------------------------------------------------------------  
     ij = 0  
     DO j=1, nb(2)  
        DO i=1, nb(1)  
           ij = ij+1  
           y(ij) = EXP(x(i, j))  
        ENDDO  
     ENDDO  
   
     nbo = nb(1)*nb(2)  
     ma_exp_r21 = 0  
     !----------------------  
   END FUNCTION ma_exp_r21  
   
   !************************************************  
   
   INTEGER FUNCTION ma_alog_r21(nb, x, nbo, y)  
     INTEGER :: nb(2), nbo, i, j, ij  
     REAL :: x(nb(1), nb(2)), y(nbo)  
     !---------------------------------------------------------------------  
     ij = 0  
     DO j=1, nb(2)  
        DO i=1, nb(1)  
           ij = ij+1  
           y(ij) = ALOG(x(i, j))  
        ENDDO  
     ENDDO  
   
     nbo = nb(1)*nb(2)  
     ma_alog_r21 = 0  
     !-----------------------  
   END FUNCTION ma_alog_r21  
   
   !************************************************  
   
   INTEGER FUNCTION ma_sqrt_r21(nb, x, nbo, y)  
     INTEGER :: nb(2), nbo, i, j, ij  
     REAL :: x(nb(1), nb(2)), y(nbo)  
     !---------------------------------------------------------------------  
     ij = 0  
     DO j=1, nb(2)  
        DO i=1, nb(1)  
           ij = ij+1  
           y(ij) = SQRT(x(i, j))  
        ENDDO  
     ENDDO  
   
     nbo = nb(1)*nb(2)  
     ma_sqrt_r21 = 0  
     !-----------------------  
   END FUNCTION ma_sqrt_r21  
   
   !************************************************  
   
   INTEGER FUNCTION ma_abs_r21(nb, x, nbo, y)  
     INTEGER :: nb(2), nbo, i, j, ij  
     REAL :: x(nb(1), nb(2)), y(nbo)  
     !---------------------------------------------------------------------  
     ij = 0  
     DO j=1, nb(2)  
        DO i=1, nb(1)  
           ij = ij+1  
           y(ij) = ABS(x(i, j))  
        ENDDO  
     ENDDO  
   
     nbo = nb(1)*nb(2)  
     ma_abs_r21 = 0  
     !----------------------  
   END FUNCTION ma_abs_r21  
   
   !************************************************  
   
   INTEGER FUNCTION ma_chs_r21(nb, x, nbo, y)  
     INTEGER :: nb(2), nbo, i, j, ij  
     REAL :: x(nb(1), nb(2)), y(nbo)  
     !---------------------------------------------------------------------  
     ij = 0  
     DO j=1, nb(2)  
        DO i=1, nb(1)  
           ij = ij+1  
           y(ij) = x(i, j)*(-1.)  
        ENDDO  
     ENDDO  
   
     nbo = nb(1)*nb(2)  
     ma_chs_r21 = 0  
     !----------------------  
   END FUNCTION ma_chs_r21  
   
   !************************************************  
   
   INTEGER FUNCTION ma_cels_r21(nb, x, nbo, y)  
     INTEGER :: nb(2), nbo, i, j, ij  
     REAL :: x(nb(1), nb(2)), y(nbo)  
     !---------------------------------------------------------------------  
     ij = 0  
     DO j=1, nb(2)  
        DO i=1, nb(1)  
           ij = ij+1  
           y(ij) = x(i, j)-273.15  
        ENDDO  
     ENDDO  
   
     nbo = nb(1)*nb(2)  
     ma_cels_r21 = 0  
     !-----------------------  
   END FUNCTION ma_cels_r21  
   
   !************************************************  
   
   INTEGER FUNCTION ma_kelv_r21(nb, x, nbo, y)  
     INTEGER :: nb(2), nbo, i, j, ij  
     REAL :: x(nb(1), nb(2)), y(nbo)  
     !---------------------------------------------------------------------  
     ij = 0  
     DO j=1, nb(2)  
        DO i=1, nb(1)  
           ij = ij+1  
           y(ij) = x(i, j)+273.15  
        ENDDO  
     ENDDO  
   
     nbo = nb(1)*nb(2)  
     ma_kelv_r21 = 0  
326      !-----------------------      !-----------------------
327    END FUNCTION ma_kelv_r21    END SUBROUTINE mathop_r21
328      !
329    !************************************************    SUBROUTINE mathop_r31(fun, nb, work_in, miss_val, nb_index, nindex, scal, &
330           nb_max, work_out)
331    INTEGER FUNCTION ma_deg_r21(nb, x, nbo, y)  
332      INTEGER :: nb(2), nbo, i, j, ij      ! This subroutines gives an interface to the various operations
333      REAL :: x(nb(1), nb(2)), y(nbo)      ! which are allowed. The interface is general enough to allow its use
334      !---------------------------------------------------------------------      ! for other cases.
335      ij = 0  
336      DO j=1, nb(2)      ! INPUT
337         DO i=1, nb(1)  
338            ij = ij+1      ! fun      : function to be applied to the vector of data
339            y(ij) = x(i, j)*57.29577951      ! nb       : Length of input vector
340         ENDDO      ! work_in  : Input vector of data (REAL)
341      ENDDO      ! miss_val : The value of the missing data flag (it has to be a
342        !            maximum value, in f90 : huge( a real ))
343      nbo = nb(1)*nb(2)      ! nb_index : Length of index vector
344      ma_deg_r21 = 0      ! nindex   : Vector of indices
345      !----------------------      ! scal     : A scalar value for vector/scalar operations
346    END FUNCTION ma_deg_r21      ! nb_max   : maximum length of output vector
   
   !************************************************  
   
   INTEGER FUNCTION ma_rad_r21(nb, x, nbo, y)  
     INTEGER :: nb(2), nbo, i, j, ij  
     REAL :: x(nb(1), nb(2)), y(nbo)  
     !---------------------------------------------------------------------  
     ij = 0  
     DO j=1, nb(2)  
        DO i=1, nb(1)  
           ij = ij+1  
           y(ij) = x(i, j)*0.01745329252  
        ENDDO  
     ENDDO  
   
     nbo = nb(1)*nb(2)  
     ma_rad_r21 = 0  
     !----------------------  
   END FUNCTION ma_rad_r21  
   
   !************************************************  
   
   INTEGER FUNCTION ma_ident_r21(nb, x, nbo, y)  
     INTEGER :: nb(2), nbo, i, j, ij  
     REAL :: x(nb(1), nb(2)), y(nbo)  
     !---------------------------------------------------------------------  
     ij = 0  
     DO j=1, nb(2)  
        DO i=1, nb(1)  
           ij = ij+1  
           y(ij) = x(i, j)  
        ENDDO  
     ENDDO  
   
     nbo = nb(1)*nb(2)  
     ma_ident_r21 = 0  
     !------------------------  
   END FUNCTION ma_ident_r21  
   !-  
   !=== OPERATIONS (two argument)  
   !-  
   INTEGER FUNCTION ma_add_r21(nb, x, s, nbo, y)  
     INTEGER :: nb(2), nbo  
     REAL :: x(nb(1), nb(2)), s, y(nbo)  
   
     INTEGER :: i, j, ij  
     !---------------------------------------------------------------------  
     ij = 0  
     DO j=1, nb(2)  
        DO i=1, nb(1)  
           ij = ij+1  
           y(ij) = x(i, j)+s  
        ENDDO  
     ENDDO  
   
     nbo = nb(1)*nb(2)  
     ma_add_r21 = 0  
     !----------------------  
   END FUNCTION ma_add_r21  
   
   !************************************************  
   
   INTEGER FUNCTION ma_sub_r21(nb, x, s, nbo, y)  
     INTEGER :: nb(2), nbo  
     REAL :: x(nb(1), nb(2)), s, y(nbo)  
   
     INTEGER :: i, j, ij  
     !---------------------------------------------------------------------  
     ij = 0  
     DO j=1, nb(2)  
        DO i=1, nb(1)  
           ij = ij+1  
           y(ij) = x(i, j)-s  
        ENDDO  
     ENDDO  
   
     nbo = nb(1)*nb(2)  
     ma_sub_r21 = 0  
     !----------------------  
   END FUNCTION ma_sub_r21  
   
   !************************************************  
   
   INTEGER FUNCTION ma_subi_r21(nb, x, s, nbo, y)  
     INTEGER :: nb(2), nbo  
     REAL :: x(nb(1), nb(2)), s, y(nbo)  
   
     INTEGER :: i, j, ij  
     !---------------------------------------------------------------------  
     ij = 0  
     DO j=1, nb(2)  
        DO i=1, nb(1)  
           ij = ij+1  
           y(ij) =  s-x(i, j)  
        ENDDO  
     ENDDO  
   
     nbo = nb(1)*nb(2)  
     ma_subi_r21 = 0  
     !-----------------------  
   END FUNCTION ma_subi_r21  
   
   !************************************************  
   
   INTEGER FUNCTION ma_mult_r21(nb, x, s, nbo, y)  
     INTEGER :: nb(2), nbo  
     REAL :: x(nb(1), nb(2)), s, y(nbo)  
   
     INTEGER :: i, j, ij  
     !---------------------------------------------------------------------  
     ij = 0  
     DO j=1, nb(2)  
        DO i=1, nb(1)  
           ij = ij+1  
           y(ij) = x(i, j)*s  
        ENDDO  
     ENDDO  
   
     nbo = nb(1)*nb(2)  
     ma_mult_r21 = 0  
     !-----------------------  
   END FUNCTION ma_mult_r21  
   
   !************************************************  
   
   INTEGER FUNCTION ma_div_r21(nb, x, s, nbo, y)  
     INTEGER :: nb(2), nbo  
     REAL :: x(nb(1), nb(2)), s, y(nbo)  
   
     INTEGER :: i, j, ij  
     !---------------------------------------------------------------------  
     ij = 0  
     DO j=1, nb(2)  
        DO i=1, nb(1)  
           ij = ij+1  
           y(ij) = x(i, j)/s  
        ENDDO  
     ENDDO  
   
     nbo = nb(1)*nb(2)  
     ma_div_r21 = 0  
     !----------------------  
   END FUNCTION ma_div_r21  
   
   !************************************************  
   
   INTEGER FUNCTION ma_divi_r21(nb, x, s, nbo, y)  
     INTEGER :: nb(2), nbo  
     REAL :: x(nb(1), nb(2)), s, y(nbo)  
   
     INTEGER :: i, j, ij  
     !---------------------------------------------------------------------  
     ij = 0  
     DO j=1, nb(2)  
        DO i=1, nb(1)  
           ij = ij+1  
           y(ij) = s/x(i, j)  
        ENDDO  
     ENDDO  
   
     nbo = nb(1)*nb(2)  
     ma_divi_r21 = 0  
     !-----------------------  
   END FUNCTION ma_divi_r21  
   
   !************************************************  
   
   INTEGER FUNCTION ma_power_r21(nb, x, s, nbo, y)  
     INTEGER :: nb(2), nbo  
     REAL :: x(nb(1), nb(2)), s, y(nbo)  
   
     INTEGER :: i, j, ij  
     !---------------------------------------------------------------------  
     ij = 0  
     DO j=1, nb(2)  
        DO i=1, nb(1)  
           ij = ij+1  
           y(ij) = x(i, j) ** s  
        ENDDO  
     ENDDO  
   
     nbo = nb(1)*nb(2)  
     ma_power_r21 = 0  
     !------------------------  
   END FUNCTION ma_power_r21  
   
   !************************************************  
   
   INTEGER FUNCTION ma_fumin_r21(nb, x, s, nbo, y)  
     INTEGER :: nb(2), nbo  
     REAL :: x(nb(1), nb(2)), s, y(nbo)  
   
     INTEGER :: i, j, ij  
     !---------------------------------------------------------------------  
     ij = 0  
     DO j=1, nb(2)  
        DO i=1, nb(1)  
           ij = ij+1  
           y(ij) = MIN(x(i, j), s)  
        ENDDO  
     ENDDO  
   
     nbo = nb(1)*nb(2)  
     ma_fumin_r21 = 0  
     !------------------------  
   END FUNCTION ma_fumin_r21  
   
   !************************************************  
   
   INTEGER FUNCTION ma_fumax_r21(nb, x, s, nbo, y)  
     INTEGER :: nb(2), nbo  
     REAL :: x(nb(1), nb(2)), s, y(nbo)  
   
     INTEGER :: i, j, ij  
     !---------------------------------------------------------------------  
     ij = 0  
     DO j=1, nb(2)  
        DO i=1, nb(1)  
           ij = ij+1  
           y(ij) = MAX(x(i, j), s)  
        ENDDO  
     ENDDO  
   
     nbo = nb(1)*nb(2)  
     ma_fumax_r21 = 0  
     !------------------------  
   END FUNCTION ma_fumax_r21  
   
   !************************************************  
   
   INTEGER FUNCTION ma_fuscat_r21(nb, x, nbi, ind, miss_val, nbo, y)  
     INTEGER :: nb(2), nbo, nbi  
     INTEGER :: ind(nbi)  
     REAL :: x(nb(1), nb(2)), miss_val, y(nbo)  
   
     INTEGER :: i, j, ij, ii, ipos  
     !---------------------------------------------------------------------  
     ma_fuscat_r21 = 0  
   
     y(1:nbo) = miss_val  
   
     IF (nbi <= nb(1)*nb(2)) THEN  
        ipos = 0  
        DO ij=1, nbi  
           IF (ind(ij) <= nbo .AND. ind(ij) > 0) THEN  
              ipos = ipos+1  
              j = ((ipos-1)/nb(1))+1  
              i = (ipos-(j-1)*nb(1))  
              y(ind(ij)) = x(i, j)  
           ELSE  
              IF (ind(ij) > nbo) ma_fuscat_r21  = ma_fuscat_r21+1  
           ENDIF  
        ENDDO  
        !-- Repeat the data if needed  
        IF (MINVAL(ind) < 0) THEN  
           DO i=1, nbi  
              IF (ind(i) <= 0) THEN  
                 DO ii=1, ABS(ind(i))-1  
                    IF (ind(i+1)+ii <= nbo) THEN  
                       y(ind(i+1)+ii) = y(ind(i+1))  
                    ELSE  
                       ma_fuscat_r21  = ma_fuscat_r21+1  
                    ENDIF  
                 ENDDO  
              ENDIF  
           ENDDO  
        ENDIF  
     ELSE  
        ma_fuscat_r21  = 1  
     ENDIF  
     !-------------------------  
   END FUNCTION ma_fuscat_r21  
   
   !************************************************  
   
   INTEGER FUNCTION ma_fugath_r21(nb, x, nbi, ind, miss_val, nbo, y)  
     INTEGER :: nb(2), nbo, nbi  
     INTEGER :: ind(nbi)  
     REAL :: x(nb(1), nb(2)), miss_val, y(nbo)  
   
     INTEGER :: i, j, ij, ipos  
     !---------------------------------------------------------------------  
     IF (nbi <= nbo) THEN  
        ma_fugath_r21 = 0  
        y(1:nbo) = miss_val  
        ipos = 0  
        DO ij=1, nbi  
           IF (ipos+1 <= nbo) THEN  
              IF (ind(ij) > 0) THEN  
                 j = ((ind(ij)-1)/nb(1))+1  
                 i = (ind(ij)-(j-1)*nb(1))  
                 ipos = ipos+1  
                 y(ipos) = x(i, j)  
              ENDIF  
           ELSE  
              IF (ipos+1 > nbo) ma_fugath_r21  = ma_fugath_r21+1  
           ENDIF  
        ENDDO  
     ELSE  
        ma_fugath_r21 = 1  
     ENDIF  
     nbo = ipos  
     !-------------------------  
   END FUNCTION ma_fugath_r21  
   
   !************************************************  
   
   INTEGER FUNCTION ma_fufill_r21(nb, x, nbi, ind, miss_val, nbo, y)  
     INTEGER :: nb(2), nbo, nbi  
     INTEGER :: ind(nbi)  
     REAL :: x(nb(1), nb(2)), miss_val, y(nbo)  
   
     INTEGER :: i, j, ij, ii, ipos  
     !---------------------------------------------------------------------  
     ma_fufill_r21 = 0  
   
     IF (nbi <= nb(1)*nb(2)) THEN  
        ipos = 0  
        DO ij=1, nbi  
           IF (ind(ij) <= nbo .AND. ind(ij) > 0) THEN  
              ipos = ipos+1  
              j = ((ipos-1)/nb(1))+1  
              i = (ipos-(j-1)*nb(1))  
              y(ind(ij)) = x(i, j)  
           ELSE  
              IF (ind(ij) > nbo) ma_fufill_r21  = ma_fufill_r21+1  
           ENDIF  
        ENDDO  
        !-- Repeat the data if needed  
        IF (MINVAL(ind) < 0) THEN  
           DO i=1, nbi  
              IF (ind(i) <= 0) THEN  
                 DO ii=1, ABS(ind(i))-1  
                    IF (ind(i+1)+ii <= nbo) THEN  
                       y(ind(i+1)+ii) = y(ind(i+1))  
                    ELSE  
                       ma_fufill_r21  = ma_fufill_r21+1  
                    ENDIF  
                 ENDDO  
              ENDIF  
           ENDDO  
        ENDIF  
     ELSE  
        ma_fufill_r21  = 1  
     ENDIF  
     !-------------------------  
   END FUNCTION ma_fufill_r21  
   
   !************************************************  
   
   INTEGER FUNCTION ma_fucoll_r21(nb, x, nbi, ind, miss_val, nbo, y)  
     INTEGER :: nb(2), nbo, nbi  
     INTEGER :: ind(nbi)  
     REAL :: x(nb(1), nb(2)), miss_val, y(nbo)  
   
     INTEGER :: i, j, ij, ipos  
     !---------------------------------------------------------------------  
     IF (nbi <= nbo) THEN  
        ma_fucoll_r21 = 0  
        ipos = 0  
        DO ij=1, nbi  
           IF (ipos+1 <= nbo) THEN  
              IF (ind(ij) > 0) THEN  
                 j = ((ind(ij)-1)/nb(1))+1  
                 i = (ind(ij)-(j-1)*nb(1))  
                 ipos = ipos+1  
                 y(ipos) = x(i, j)  
              ENDIF  
           ELSE  
              IF (ipos+1 > nbo) ma_fucoll_r21  = ma_fucoll_r21+1  
           ENDIF  
        ENDDO  
     ELSE  
        ma_fucoll_r21 = 1  
     ENDIF  
     nbo = ipos  
     !-------------------------  
   END FUNCTION ma_fucoll_r21  
   
   !************************************************  
   
   INTEGER FUNCTION ma_fuundef_r21(nb, x, nbi, ind, miss_val, nbo, y)  
     INTEGER :: nb(2), nbo, nbi  
     INTEGER :: ind(nbi)  
     REAL :: x(nb(1), nb(2)), miss_val, y(nbo)  
   
     INTEGER :: i, j, ij  
     !---------------------------------------------------------------------  
     IF (nbi <= nbo .AND. nbo == nb(1)*nb(2)) THEN  
        ma_fuundef_r21 = 0  
        DO ij=1, nbo  
           j = ((ij-1)/nb(1))+1  
           i = (ij-(j-1)*nb(1))  
           y(ij) = x(i, j)  
        ENDDO  
        DO i=1, nbi  
           IF (ind(i) <= nbo .AND. ind(i) > 0) THEN  
              y(ind(i)) =  miss_val  
           ELSE  
              IF (ind(i) > nbo) ma_fuundef_r21 = ma_fuundef_r21+1  
           ENDIF  
        ENDDO  
     ELSE  
        ma_fuundef_r21 = 1  
     ENDIF  
     !--------------------------  
   END FUNCTION ma_fuundef_r21  
   
   !************************************************  
   
   INTEGER FUNCTION ma_fuonly_r21(nb, x, nbi, ind, miss_val, nbo, y)  
     INTEGER :: nb(2), nbo, nbi  
     INTEGER :: ind(nbi)  
     REAL :: x(nb(1), nb(2)), miss_val, y(nbo)  
   
     INTEGER :: i, j, ij  
     !---------------------------------------------------------------------  
     IF (     (nbi <= nbo).AND.(nbo == nb(1)*nb(2)) &  
          &    .AND.ALL(ind(1:nbi) <= nbo) ) THEN  
        ma_fuonly_r21 = 0  
        y(1:nbo) = miss_val  
        DO ij=1, nbi  
           IF (ind(ij) > 0) THEN  
              j = ((ind(ij)-1)/nb(1))+1  
              i = (ind(ij)-(j-1)*nb(1))  
              y(ind(ij)) =  x(i, j)  
           ENDIF  
        ENDDO  
     ELSE  
        ma_fuonly_r21 = 1  
     ENDIF  
     !-------------------------  
   END FUNCTION ma_fuonly_r21  
   
   !************************************************  
   
   !************************************************  
   
   SUBROUTINE mathop_r31 &  
        &  (fun, nb, work_in, miss_val, nb_index, nindex, scal, nb_max, work_out)  
     !- This subroutines gives an interface to the various operations  
     !- which are allowed. The interface is general enough to allow its use  
     !- for other cases.  
   
     !- INPUT  
   
     !- fun      : function to be applied to the vector of data  
     !- nb       : Length of input vector  
     !- work_in  : Input vector of data (REAL)  
     !- miss_val : The value of the missing data flag (it has to be a  
     !-            maximum value, in f90 : huge( a real ))  
     !- nb_index : Length of index vector  
     !- nindex   : Vector of indices  
     !- scal     : A scalar value for vector/scalar operations  
     !- nb_max   : maximum length of output vector  
347    
348      !- OUTPUT      ! OUTPUT
349    
350      !- nb_max   : Actual length of output variable      ! nb_max   : Actual length of output variable
351      !- work_out : Output vector after the operation was applied      ! work_out : Output vector after the operation was applied
352      USE errioipsl, ONLY : histerr      USE errioipsl, ONLY : histerr
353        USE mathop2, ONLY : ma_abs_r31, ma_acos_r31, ma_add_r31, ma_alog_r31, &
354             ma_asin_r31, ma_atan_r31, ma_cels_r31, ma_chs_r31, ma_cos_r31, &
355             ma_deg_r31, ma_divi_r31, ma_div_r31, ma_exp_r31, ma_fucoll_r31, &
356             ma_fufill_r31, ma_fugath_r31, ma_fumax_r31, ma_fumin_r31, &
357             ma_fuonly_r31, ma_fuscat_r31, ma_fuundef_r31, ma_ident_r31, &
358             ma_kelv_r31, ma_mult_r31, ma_power_r31, ma_rad_r31, ma_sin_r31, &
359             ma_sqrt_r31, ma_subi_r31, ma_sub_r31, ma_tan_r31
360    
361      CHARACTER(LEN=7) :: fun      CHARACTER(LEN=7) :: fun
362      INTEGER :: nb(3), nb_max, nb_index      INTEGER :: nb(3), nb_max, nb_index
# Line 1640  CONTAINS Line 367  CONTAINS
367      INTEGER :: ierr      INTEGER :: ierr
368    
369      INTRINSIC SIN, COS, TAN, ASIN, ACOS, ATAN, EXP, ALOG, SQRT, ABS      INTRINSIC SIN, COS, TAN, ASIN, ACOS, ATAN, EXP, ALOG, SQRT, ABS
370      !---------------------------------------------------------------------      !--------------------------------------------------------------------
371      ierr = 0      ierr = 0
372    
373      IF (scal >= miss_val-1.) THEN      IF (scal >= miss_val-1.) THEN
# Line 1680  CONTAINS Line 407  CONTAINS
407               ierr = ma_ident_r31(nb, work_in, nb_max, work_out)               ierr = ma_ident_r31(nb, work_in, nb_max, work_out)
408            CASE DEFAULT            CASE DEFAULT
409               CALL histerr(3, "mathop", &               CALL histerr(3, "mathop", &
410                    &        'scalar variable undefined and no indexing', &                    'scalar variable undefined and no indexing', &
411                    &        'but still unknown function', fun)                    'but still unknown function', fun)
412            END SELECT            END SELECT
413            IF (ierr > 0) THEN            IF (ierr > 0) THEN
414               CALL histerr(3, "mathop", &               CALL histerr(3, "mathop", &
415                    &        'Error while executing a simple function', fun, ' ')                    'Error while executing a simple function', fun, ' ')
416            ENDIF            ENDIF
417         ELSE         ELSE
418            SELECT CASE (fun)            SELECT CASE (fun)
419            CASE('gather')            CASE('gather')
420               ierr = ma_fugath_r31(nb, work_in, nb_index, nindex, &               ierr = ma_fugath_r31(nb, work_in, nb_index, nindex, &
421                    &                           miss_val, nb_max, work_out)                    miss_val, nb_max, work_out)
422            CASE('scatter')            CASE('scatter')
423               IF (nb_index > (nb(1)*nb(2)*nb(3))) THEN               IF (nb_index > (nb(1)*nb(2)*nb(3))) THEN
424                  work_out(1:nb_max) = miss_val                  work_out(1:nb_max) = miss_val
425                  ierr=1                  ierr=1
426               ELSE               ELSE
427                  ierr = ma_fuscat_r31(nb, work_in, nb_index, nindex, &                  ierr = ma_fuscat_r31(nb, work_in, nb_index, nindex, &
428                       &                             miss_val, nb_max, work_out)                       miss_val, nb_max, work_out)
429               ENDIF               ENDIF
430            CASE('coll')            CASE('coll')
431               ierr = ma_fucoll_r31(nb, work_in, nb_index, nindex, &               ierr = ma_fucoll_r31(nb, work_in, nb_index, nindex, &
432                    &                           miss_val, nb_max, work_out)                    miss_val, nb_max, work_out)
433            CASE('fill')            CASE('fill')
434               ierr = ma_fufill_r31(nb, work_in, nb_index, nindex, &               ierr = ma_fufill_r31(nb, work_in, nb_index, nindex, &
435                    &                           miss_val, nb_max, work_out)                    miss_val, nb_max, work_out)
436            CASE('undef')            CASE('undef')
437               ierr = ma_fuundef_r31(nb, work_in, nb_index, nindex, &               ierr = ma_fuundef_r31(nb, work_in, nb_index, nindex, &
438                    &                           miss_val, nb_max, work_out)                    miss_val, nb_max, work_out)
439            CASE('only')            CASE('only')
440               ierr = ma_fuonly_r31(nb, work_in, nb_index, nindex, &               ierr = ma_fuonly_r31(nb, work_in, nb_index, nindex, &
441                    &                           miss_val, nb_max, work_out)                    miss_val, nb_max, work_out)
442            CASE DEFAULT            CASE DEFAULT
443               CALL histerr(3, "mathop", &               CALL histerr(3, "mathop", &
444                    &        'scalar variable undefined and indexing', &                    'scalar variable undefined and indexing', &
445                    &        'was requested but with unknown function', fun)                    'was requested but with unknown function', fun)
446            END SELECT            END SELECT
447            IF (ierr > 0) THEN            IF (ierr > 0) THEN
448               CALL histerr(3, "mathop_r31", &               CALL histerr(3, "mathop_r31", &
449                    &        'Error while executing an indexing function', fun, ' ')                    'Error while executing an indexing function', fun, ' ')
450            ENDIF            ENDIF
451         ENDIF         ENDIF
452      ELSE      ELSE
# Line 1744  CONTAINS Line 471  CONTAINS
471            ierr = ma_power_r31(nb, work_in, scal, nb_max, work_out)            ierr = ma_power_r31(nb, work_in, scal, nb_max, work_out)
472         CASE DEFAULT         CASE DEFAULT
473            CALL histerr(3, "mathop", &            CALL histerr(3, "mathop", &
474                 &      'Unknown operation with a scalar', fun, ' ')                 'Unknown operation with a scalar', fun, ' ')
475         END SELECT         END SELECT
476         IF (ierr > 0) THEN         IF (ierr > 0) THEN
477            CALL histerr(3, "mathop", &            CALL histerr(3, "mathop", &
478                 &      'Error while executing a scalar function', fun, ' ')                 'Error while executing a scalar function', fun, ' ')
479         ENDIF         ENDIF
480      ENDIF      ENDIF
     !------------------------  
   END SUBROUTINE mathop_r31  
   !-  
   !=== FUNCTIONS (only one argument)  
   !-  
   INTEGER FUNCTION ma_sin_r31(nb, x, nbo, y)  
     INTEGER :: nb(3), nbo, i, j, k, ij  
     REAL :: x(nb(1), nb(2), nb(3)), y(nbo)  
     !---------------------------------------------------------------------  
     ij = 0  
     DO k=1, nb(3)  
        DO j=1, nb(2)  
           DO i=1, nb(1)  
              ij = ij+1  
              y(ij) = SIN(x(i, j, k))  
           ENDDO  
        ENDDO  
     ENDDO  
   
     nbo = nb(1)*nb(2)*nb(3)  
     ma_sin_r31 = 0  
     !----------------------  
   END FUNCTION ma_sin_r31  
   
   !************************************************  
   
   INTEGER FUNCTION ma_cos_r31(nb, x, nbo, y)  
     INTEGER :: nb(3), nbo, i, j, k, ij  
     REAL :: x(nb(1), nb(2), nb(3)), y(nbo)  
     !---------------------------------------------------------------------  
     ij = 0  
     DO k=1, nb(3)  
        DO j=1, nb(2)  
           DO i=1, nb(1)  
              ij = ij+1  
              y(ij) = COS(x(i, j, k))  
           ENDDO  
        ENDDO  
     ENDDO  
   
     nbo = nb(1)*nb(2)*nb(3)  
     ma_cos_r31 = 0  
     !----------------------  
   END FUNCTION ma_cos_r31  
   
   !************************************************  
   
   INTEGER FUNCTION ma_tan_r31(nb, x, nbo, y)  
     INTEGER :: nb(3), nbo, i, j, k, ij  
     REAL :: x(nb(1), nb(2), nb(3)), y(nbo)  
     !---------------------------------------------------------------------  
     ij = 0  
     DO k=1, nb(3)  
        DO j=1, nb(2)  
           DO i=1, nb(1)  
              ij = ij+1  
              y(ij) = TAN(x(i, j, k))  
           ENDDO  
        ENDDO  
     ENDDO  
   
     nbo = nb(1)*nb(2)*nb(3)  
     ma_tan_r31 = 0  
     !----------------------  
   END FUNCTION ma_tan_r31  
   
   !************************************************  
   
   INTEGER FUNCTION ma_asin_r31(nb, x, nbo, y)  
     INTEGER :: nb(3), nbo, i, j, k, ij  
     REAL :: x(nb(1), nb(2), nb(3)), y(nbo)  
     !---------------------------------------------------------------------  
     ij = 0  
     DO k=1, nb(3)  
        DO j=1, nb(2)  
           DO i=1, nb(1)  
              ij = ij+1  
              y(ij) = ASIN(x(i, j, k))  
           ENDDO  
        ENDDO  
     ENDDO  
   
     nbo = nb(1)*nb(2)*nb(3)  
     ma_asin_r31 = 0  
     !-----------------------  
   END FUNCTION ma_asin_r31  
   
   !************************************************  
   
   INTEGER FUNCTION ma_acos_r31(nb, x, nbo, y)  
     INTEGER :: nb(3), nbo, i, j, k, ij  
     REAL :: x(nb(1), nb(2), nb(3)), y(nbo)  
     !---------------------------------------------------------------------  
     ij = 0  
     DO k=1, nb(3)  
        DO j=1, nb(2)  
           DO i=1, nb(1)  
              ij = ij+1  
              y(ij) = ACOS(x(i, j, k))  
           ENDDO  
        ENDDO  
     ENDDO  
   
     nbo = nb(1)*nb(2)*nb(3)  
     ma_acos_r31 = 0  
     !-----------------------  
   END FUNCTION ma_acos_r31  
   
   !************************************************  
   
   INTEGER FUNCTION ma_atan_r31(nb, x, nbo, y)  
     INTEGER :: nb(3), nbo, i, j, k, ij  
     REAL :: x(nb(1), nb(2), nb(3)), y(nbo)  
     !---------------------------------------------------------------------  
     ij = 0  
     DO k=1, nb(3)  
        DO j=1, nb(2)  
           DO i=1, nb(1)  
              ij = ij+1  
              y(ij) = ATAN(x(i, j, k))  
           ENDDO  
        ENDDO  
     ENDDO  
   
     nbo = nb(1)*nb(2)*nb(3)  
     ma_atan_r31 = 0  
     !-----------------------  
   END FUNCTION ma_atan_r31  
   
   !************************************************  
481    
482    INTEGER FUNCTION ma_exp_r31(nb, x, nbo, y)    END SUBROUTINE mathop_r31
     INTEGER :: nb(3), nbo, i, j, k, ij  
     REAL :: x(nb(1), nb(2), nb(3)), y(nbo)  
     !---------------------------------------------------------------------  
     ij = 0  
     DO k=1, nb(3)  
        DO j=1, nb(2)  
           DO i=1, nb(1)  
              ij = ij+1  
              y(ij) = EXP(x(i, j, k))  
           ENDDO  
        ENDDO  
     ENDDO  
   
     nbo = nb(1)*nb(2)*nb(3)  
     ma_exp_r31 = 0  
     !----------------------  
   END FUNCTION ma_exp_r31  
   
   !************************************************  
   
   INTEGER FUNCTION ma_alog_r31(nb, x, nbo, y)  
     INTEGER :: nb(3), nbo, i, j, k, ij  
     REAL :: x(nb(1), nb(2), nb(3)), y(nbo)  
     !---------------------------------------------------------------------  
     ij = 0  
     DO k=1, nb(3)  
        DO j=1, nb(2)  
           DO i=1, nb(1)  
              ij = ij+1  
              y(ij) = ALOG(x(i, j, k))  
           ENDDO  
        ENDDO  
     ENDDO  
   
     nbo = nb(1)*nb(2)*nb(3)  
     ma_alog_r31 = 0  
     !-----------------------  
   END FUNCTION ma_alog_r31  
   
   !************************************************  
   
   INTEGER FUNCTION ma_sqrt_r31(nb, x, nbo, y)  
     INTEGER :: nb(3), nbo, i, j, k, ij  
     REAL :: x(nb(1), nb(2), nb(3)), y(nbo)  
     !---------------------------------------------------------------------  
     ij = 0  
     DO k=1, nb(3)  
        DO j=1, nb(2)  
           DO i=1, nb(1)  
              ij = ij+1  
              y(ij) = SQRT(x(i, j, k))  
           ENDDO  
        ENDDO  
     ENDDO  
   
     nbo = nb(1)*nb(2)*nb(3)  
     ma_sqrt_r31 = 0  
     !-----------------------  
   END FUNCTION ma_sqrt_r31  
   
   !************************************************  
   
   INTEGER FUNCTION ma_abs_r31(nb, x, nbo, y)  
     INTEGER :: nb(3), nbo, i, j, k, ij  
     REAL :: x(nb(1), nb(2), nb(3)), y(nbo)  
     !---------------------------------------------------------------------  
     ij = 0  
     DO k=1, nb(3)  
        DO j=1, nb(2)  
           DO i=1, nb(1)  
              ij = ij+1  
              y(ij) = ABS(x(i, j, k))  
           ENDDO  
        ENDDO  
     ENDDO  
   
     nbo = nb(1)*nb(2)*nb(3)  
     ma_abs_r31 = 0  
     !----------------------  
   END FUNCTION ma_abs_r31  
   
   !************************************************  
   
   INTEGER FUNCTION ma_chs_r31(nb, x, nbo, y)  
     INTEGER :: nb(3), nbo, i, j, k, ij  
     REAL :: x(nb(1), nb(2), nb(3)), y(nbo)  
     !---------------------------------------------------------------------  
     ij = 0  
     DO k=1, nb(3)  
        DO j=1, nb(2)  
           DO i=1, nb(1)  
              ij = ij+1  
              y(ij) = x(i, j, k)*(-1.)  
           ENDDO  
        ENDDO  
     ENDDO  
   
     nbo = nb(1)*nb(2)*nb(3)  
     ma_chs_r31 = 0  
     !----------------------  
   END FUNCTION ma_chs_r31  
   
   !************************************************  
   
   INTEGER FUNCTION ma_cels_r31(nb, x, nbo, y)  
     INTEGER :: nb(3), nbo, i, j, k, ij  
     REAL :: x(nb(1), nb(2), nb(3)), y(nbo)  
     !---------------------------------------------------------------------  
     ij = 0  
     DO k=1, nb(3)  
        DO j=1, nb(2)  
           DO i=1, nb(1)  
              ij = ij+1  
              y(ij) = x(i, j, k)-273.15  
           ENDDO  
        ENDDO  
     ENDDO  
   
     nbo = nb(1)*nb(2)*nb(3)  
     ma_cels_r31 = 0  
     !-----------------------  
   END FUNCTION ma_cels_r31  
   
   !************************************************  
   
   INTEGER FUNCTION ma_kelv_r31(nb, x, nbo, y)  
     INTEGER :: nb(3), nbo, i, j, k, ij  
     REAL :: x(nb(1), nb(2), nb(3)), y(nbo)  
     !---------------------------------------------------------------------  
     ij = 0  
     DO k=1, nb(3)  
        DO j=1, nb(2)  
           DO i=1, nb(1)  
              ij = ij+1  
              y(ij) = x(i, j, k)+273.15  
           ENDDO  
        ENDDO  
     ENDDO  
   
     nbo = nb(1)*nb(2)*nb(3)  
     ma_kelv_r31 = 0  
     !-----------------------  
   END FUNCTION ma_kelv_r31  
   
   !************************************************  
   
   INTEGER FUNCTION ma_deg_r31(nb, x, nbo, y)  
     INTEGER :: nb(3), nbo, i, j, k, ij  
     REAL :: x(nb(1), nb(2), nb(3)), y(nbo)  
     !---------------------------------------------------------------------  
     ij = 0  
     DO k=1, nb(3)  
        DO j=1, nb(2)  
           DO i=1, nb(1)  
              ij = ij+1  
              y(ij) = x(i, j, k)*57.29577951  
           ENDDO  
        ENDDO  
     ENDDO  
   
     nbo = nb(1)*nb(2)*nb(3)  
     ma_deg_r31 = 0  
     !----------------------  
   END FUNCTION ma_deg_r31  
   
   !************************************************  
   
   INTEGER FUNCTION ma_rad_r31(nb, x, nbo, y)  
     INTEGER :: nb(3), nbo, i, j, k, ij  
     REAL :: x(nb(1), nb(2), nb(3)), y(nbo)  
     !---------------------------------------------------------------------  
     ij = 0  
     DO k=1, nb(3)  
        DO j=1, nb(2)  
           DO i=1, nb(1)  
              ij = ij+1  
              y(ij) = x(i, j, k)*0.01745329252  
           ENDDO  
        ENDDO  
     ENDDO  
   
     nbo = nb(1)*nb(2)*nb(3)  
     ma_rad_r31 = 0  
     !----------------------  
   END FUNCTION ma_rad_r31  
   
   !************************************************  
   
   INTEGER FUNCTION ma_ident_r31(nb, x, nbo, y)  
     INTEGER :: nb(3), nbo, i, j, k, ij  
     REAL :: x(nb(1), nb(2), nb(3)), y(nbo)  
     !---------------------------------------------------------------------  
     ij = 0  
     DO k=1, nb(3)  
        DO j=1, nb(2)  
           DO i=1, nb(1)  
              ij = ij+1  
              y(ij) = x(i, j, k)  
           ENDDO  
        ENDDO  
     ENDDO  
   
     nbo = nb(1)*nb(2)*nb(3)  
     ma_ident_r31 = 0  
     !------------------------  
   END FUNCTION ma_ident_r31  
   !-  
   !=== OPERATIONS (two argument)  
   !-  
   INTEGER FUNCTION ma_add_r31(nb, x, s, nbo, y)  
     INTEGER :: nb(3), nbo  
     REAL :: x(nb(1), nb(2), nb(3)), s, y(nbo)  
   
     INTEGER :: i, j, k, ij  
     !---------------------------------------------------------------------  
     ij = 0  
     DO k=1, nb(3)  
        DO j=1, nb(2)  
           DO i=1, nb(1)  
              ij = ij+1  
              y(ij) = x(i, j, k)+s  
           ENDDO  
        ENDDO  
     ENDDO  
   
     nbo = nb(1)*nb(2)*nb(3)  
     ma_add_r31 = 0  
     !----------------------  
   END FUNCTION ma_add_r31  
   
   !************************************************  
   
   INTEGER FUNCTION ma_sub_r31(nb, x, s, nbo, y)  
     INTEGER :: nb(3), nbo  
     REAL :: x(nb(1), nb(2), nb(3)), s, y(nbo)  
   
     INTEGER :: i, j, k, ij  
     !---------------------------------------------------------------------  
     ij = 0  
     DO k=1, nb(3)  
        DO j=1, nb(2)  
           DO i=1, nb(1)  
              ij = ij+1  
              y(ij) = x(i, j, k)-s  
           ENDDO  
        ENDDO  
     ENDDO  
   
     nbo = nb(1)*nb(2)*nb(3)  
     ma_sub_r31 = 0  
     !----------------------  
   END FUNCTION ma_sub_r31  
   
   !************************************************  
   
   INTEGER FUNCTION ma_subi_r31(nb, x, s, nbo, y)  
     INTEGER :: nb(3), nbo  
     REAL :: x(nb(1), nb(2), nb(3)), s, y(nbo)  
   
     INTEGER :: i, j, k, ij  
     !---------------------------------------------------------------------  
     ij = 0  
     DO k=1, nb(3)  
        DO j=1, nb(2)  
           DO i=1, nb(1)  
              ij = ij+1  
              y(ij) =  s-x(i, j, k)  
           ENDDO  
        ENDDO  
     ENDDO  
   
     nbo = nb(1)*nb(2)*nb(3)  
     ma_subi_r31 = 0  
     !-----------------------  
   END FUNCTION ma_subi_r31  
   
   !************************************************  
   
   INTEGER FUNCTION ma_mult_r31(nb, x, s, nbo, y)  
     INTEGER :: nb(3), nbo  
     REAL :: x(nb(1), nb(2), nb(3)), s, y(nbo)  
   
     INTEGER :: i, j, k, ij  
     !---------------------------------------------------------------------  
     ij = 0  
     DO k=1, nb(3)  
        DO j=1, nb(2)  
           DO i=1, nb(1)  
              ij = ij+1  
              y(ij) = x(i, j, k)*s  
           ENDDO  
        ENDDO  
     ENDDO  
   
     nbo = nb(1)*nb(2)*nb(3)  
     ma_mult_r31 = 0  
     !-----------------------  
   END FUNCTION ma_mult_r31  
   
   !************************************************  
   
   INTEGER FUNCTION ma_div_r31(nb, x, s, nbo, y)  
     INTEGER :: nb(3), nbo  
     REAL :: x(nb(1), nb(2), nb(3)), s, y(nbo)  
   
     INTEGER :: i, j, k, ij  
     !---------------------------------------------------------------------  
     ij = 0  
     DO k=1, nb(3)  
        DO j=1, nb(2)  
           DO i=1, nb(1)  
              ij = ij+1  
              y(ij) = x(i, j, k)/s  
           ENDDO  
        ENDDO  
     ENDDO  
   
     nbo = nb(1)*nb(2)*nb(3)  
     ma_div_r31 = 0  
     !----------------------  
   END FUNCTION ma_div_r31  
   
   !************************************************  
   
   INTEGER FUNCTION ma_divi_r31(nb, x, s, nbo, y)  
     INTEGER :: nb(3), nbo  
     REAL :: x(nb(1), nb(2), nb(3)), s, y(nbo)  
   
     INTEGER :: i, j, k, ij  
     !---------------------------------------------------------------------  
     ij = 0  
     DO k=1, nb(3)  
        DO j=1, nb(2)  
           DO i=1, nb(1)  
              ij = ij+1  
              y(ij) = s/x(i, j, k)  
           ENDDO  
        ENDDO  
     ENDDO  
   
     nbo = nb(1)*nb(2)*nb(3)  
     ma_divi_r31 = 0  
     !-----------------------  
   END FUNCTION ma_divi_r31  
   
   !************************************************  
   
   INTEGER FUNCTION ma_power_r31(nb, x, s, nbo, y)  
     INTEGER :: nb(3), nbo  
     REAL :: x(nb(1), nb(2), nb(3)), s, y(nbo)  
   
     INTEGER :: i, j, k, ij  
     !---------------------------------------------------------------------  
     ij = 0  
     DO k=1, nb(3)  
        DO j=1, nb(2)  
           DO i=1, nb(1)  
              ij = ij+1  
              y(ij) = x(i, j, k)**s  
           ENDDO  
        ENDDO  
     ENDDO  
   
     nbo = nb(1)*nb(2)*nb(3)  
     ma_power_r31 = 0  
     !------------------------  
   END FUNCTION ma_power_r31  
   
   !************************************************  
   
   INTEGER FUNCTION ma_fumin_r31(nb, x, s, nbo, y)  
     INTEGER :: nb(3), nbo  
     REAL :: x(nb(1), nb(2), nb(3)), s, y(nbo)  
   
     INTEGER :: i, j, k, ij  
     !---------------------------------------------------------------------  
     ij = 0  
     DO k=1, nb(3)  
        DO j=1, nb(2)  
           DO i=1, nb(1)  
              ij = ij+1  
              y(ij) = MIN(x(i, j, k), s)  
           ENDDO  
        ENDDO  
     ENDDO  
   
     nbo = nb(1)*nb(2)*nb(3)  
     ma_fumin_r31 = 0  
     !------------------------  
   END FUNCTION ma_fumin_r31  
   
   !************************************************  
   
   INTEGER FUNCTION ma_fumax_r31(nb, x, s, nbo, y)  
     INTEGER :: nb(3), nbo  
     REAL :: x(nb(1), nb(2), nb(3)), s, y(nbo)  
   
     INTEGER :: i, j, k, ij  
     !---------------------------------------------------------------------  
     ij = 0  
     DO k=1, nb(3)  
        DO j=1, nb(2)  
           DO i=1, nb(1)  
              ij = ij+1  
              y(ij) = MAX(x(i, j, k), s)  
           ENDDO  
        ENDDO  
     ENDDO  
   
     nbo = nb(1)*nb(2)*nb(3)  
     ma_fumax_r31 = 0  
     !------------------------  
   END FUNCTION ma_fumax_r31  
   
   !************************************************  
   
   INTEGER FUNCTION ma_fuscat_r31(nb, x, nbi, ind, miss_val, nbo, y)  
     INTEGER :: nb(3), nbo, nbi  
     INTEGER :: ind(nbi)  
     REAL :: x(nb(1), nb(2), nb(3)), miss_val, y(nbo)  
   
     INTEGER :: i, j, k, ij, ii, ipos, ipp, isb  
     !---------------------------------------------------------------------  
     ma_fuscat_r31 = 0  
   
     y(1:nbo) = miss_val  
   
     IF (nbi <= nb(1)*nb(2)*nb(3)) THEN  
        ipos = 0  
        isb = nb(1)*nb(2)  
        DO ij=1, nbi  
           IF (ind(ij) <= nbo .AND. ind(ij) > 0) THEN  
              ipos = ipos+1  
              k = ((ipos-1)/isb)+1  
              ipp = ipos-(k-1)*isb  
              j = ((ipp-1)/nb(1))+1  
              i = (ipp-(j-1)*nb(1))  
              y(ind(ij)) = x(i, j, k)  
           ELSE  
              IF (ind(ij) > nbo) ma_fuscat_r31  = ma_fuscat_r31+1  
           ENDIF  
        ENDDO  
        !-- Repeat the data if needed  
        IF (MINVAL(ind) < 0) THEN  
           DO i=1, nbi  
              IF (ind(i) <= 0) THEN  
                 DO ii=1, ABS(ind(i))-1  
                    IF (ind(i+1)+ii <= nbo) THEN  
                       y(ind(i+1)+ii) = y(ind(i+1))  
                    ELSE  
                       ma_fuscat_r31  = ma_fuscat_r31+1  
                    ENDIF  
                 ENDDO  
              ENDIF  
           ENDDO  
        ENDIF  
     ELSE  
        ma_fuscat_r31  = 1  
     ENDIF  
     !-------------------------  
   END FUNCTION ma_fuscat_r31  
   
   !************************************************  
   
   INTEGER FUNCTION ma_fugath_r31(nb, x, nbi, ind, miss_val, nbo, y)  
     INTEGER :: nb(3), nbo, nbi  
     INTEGER :: ind(nbi)  
     REAL :: x(nb(1), nb(2), nb(3)), miss_val, y(nbo)  
   
     INTEGER :: i, j, k, ij, ipos, ipp, isb  
     !---------------------------------------------------------------------  
     IF (nbi <= nbo) THEN  
        ma_fugath_r31 = 0  
        y(1:nbo) = miss_val  
        ipos = 0  
        isb = nb(1)*nb(2)  
        DO ij=1, nbi  
           IF (ipos+1 <= nbo) THEN  
              IF (ind(ij) > 0) THEN  
                 k = ((ind(ij)-1)/isb)+1  
                 ipp = ind(ij)-(k-1)*isb  
                 j = ((ipp-1)/nb(1))+1  
                 i = (ipp-(j-1)*nb(1))  
                 ipos = ipos+1  
                 y(ipos) = x(i, j, k)  
              ENDIF  
           ELSE  
              IF (ipos+1 > nbo) ma_fugath_r31  = ma_fugath_r31+1  
           ENDIF  
        ENDDO  
     ELSE  
        ma_fugath_r31 = 1  
     ENDIF  
     nbo = ipos  
     !-------------------------  
   END FUNCTION ma_fugath_r31  
   
   !************************************************  
   
   INTEGER FUNCTION ma_fufill_r31(nb, x, nbi, ind, miss_val, nbo, y)  
     INTEGER :: nb(3), nbo, nbi  
     INTEGER :: ind(nbi)  
     REAL :: x(nb(1), nb(2), nb(3)), miss_val, y(nbo)  
   
     INTEGER :: i, j, k, ij, ii, ipos, ipp, isb  
     !---------------------------------------------------------------------  
     ma_fufill_r31 = 0  
     IF (nbi <= nb(1)*nb(2)*nb(3)) THEN  
        ipos = 0  
        isb = nb(1)*nb(2)  
        DO ij=1, nbi  
           IF (ind(ij) <= nbo .AND. ind(ij) > 0) THEN  
              ipos = ipos+1  
              k = ((ipos-1)/isb)+1  
              ipp = ipos-(k-1)*isb  
              j = ((ipp-1)/nb(1))+1  
              i = (ipp-(j-1)*nb(1))  
              y(ind(ij)) = x(i, j, k)  
           ELSE  
              IF (ind(ij) > nbo) ma_fufill_r31  = ma_fufill_r31+1  
           ENDIF  
        ENDDO  
        !-- Repeat the data if needed  
        IF (MINVAL(ind) < 0) THEN  
           DO i=1, nbi  
              IF (ind(i) <= 0) THEN  
                 DO ii=1, ABS(ind(i))-1  
                    IF (ind(i+1)+ii <= nbo) THEN  
                       y(ind(i+1)+ii) = y(ind(i+1))  
                    ELSE  
                       ma_fufill_r31  = ma_fufill_r31+1  
                    ENDIF  
                 ENDDO  
              ENDIF  
           ENDDO  
        ENDIF  
     ELSE  
        ma_fufill_r31  = 1  
     ENDIF  
     !-------------------------  
   END FUNCTION ma_fufill_r31  
   
   !************************************************  
   
   INTEGER FUNCTION ma_fucoll_r31(nb, x, nbi, ind, miss_val, nbo, y)  
     INTEGER :: nb(3), nbo, nbi  
     INTEGER :: ind(nbi)  
     REAL :: x(nb(1), nb(2), nb(3)), miss_val, y(nbo)  
   
     INTEGER :: i, j, k, ij, ipos, ipp, isb  
     !---------------------------------------------------------------------  
     IF (nbi <= nbo) THEN  
        ma_fucoll_r31 = 0  
        ipos = 0  
        isb = nb(1)*nb(2)  
        DO ij=1, nbi  
           IF (ipos+1 <= nbo) THEN  
              IF (ind(ij) > 0) THEN  
                 k = ((ind(ij)-1)/isb)+1  
                 ipp = ind(ij)-(k-1)*isb  
                 j = ((ipp-1)/nb(1))+1  
                 i = (ipp-(j-1)*nb(1))  
                 ipos = ipos+1  
                 y(ipos) = x(i, j, k)  
              ENDIF  
           ELSE  
              IF (ipos+1 > nbo) ma_fucoll_r31  = ma_fucoll_r31+1  
           ENDIF  
        ENDDO  
     ELSE  
        ma_fucoll_r31 = 1  
     ENDIF  
     nbo = ipos  
     !-------------------------  
   END FUNCTION ma_fucoll_r31  
   
   !************************************************  
   
   INTEGER FUNCTION ma_fuundef_r31(nb, x, nbi, ind, miss_val, nbo, y)  
     INTEGER :: nb(3), nbo, nbi  
     INTEGER :: ind(nbi)  
     REAL :: x(nb(1), nb(2), nb(3)), miss_val, y(nbo)  
   
     INTEGER :: i, j, k, ij, ipp, isb  
     !---------------------------------------------------------------------  
     IF (nbi <= nbo .AND. nbo == nb(1)*nb(2)*nb(3)) THEN  
        ma_fuundef_r31 = 0  
        isb = nb(1)*nb(2)  
        DO ij=1, nbo  
           k = ((ij-1)/isb)+1  
           ipp = ij-(k-1)*isb  
           j = ((ipp-1)/nb(1))+1  
           i = (ipp-(j-1)*nb(1))  
           y(ij) = x(i, j, k)  
        ENDDO  
        DO i=1, nbi  
           IF (ind(i) <= nbo .AND. ind(i) > 0) THEN  
              y(ind(i)) =  miss_val  
           ELSE  
              IF (ind(i) > nbo) ma_fuundef_r31  = ma_fuundef_r31+1  
           ENDIF  
        ENDDO  
     ELSE  
        ma_fuundef_r31 = 1  
     ENDIF  
     !--------------------------  
   END FUNCTION ma_fuundef_r31  
   
   !************************************************  
   
   INTEGER FUNCTION ma_fuonly_r31(nb, x, nbi, ind, miss_val, nbo, y)  
     INTEGER :: nb(3), nbo, nbi  
     INTEGER :: ind(nbi)  
     REAL :: x(nb(1), nb(2), nb(3)), miss_val, y(nbo)  
   
     INTEGER :: i, j, k, ij, ipp, isb  
     !---------------------------------------------------------------------  
     IF (     (nbi <= nbo).AND.(nbo == nb(1)*nb(2)*nb(3)) &  
          &    .AND.ALL(ind(1:nbi) <= nbo) ) THEN  
        ma_fuonly_r31 = 0  
        y(1:nbo) = miss_val  
        isb = nb(1)*nb(2)  
        DO ij=1, nbi  
           IF (ind(ij) > 0) THEN  
              k = ((ind(ij)-1)/isb)+1  
              ipp = ind(ij)-(k-1)*isb  
              j = ((ipp-1)/nb(1))+1  
              i = (ipp-(j-1)*nb(1))  
              y(ind(ij)) = x(i, j, k)  
           ENDIF  
        ENDDO  
     ELSE  
        ma_fuonly_r31 = 1  
     ENDIF  
     !-------------------------  
   END FUNCTION ma_fuonly_r31  
483    
484  END MODULE mathop_m  END MODULE mathop_m
Legend:
Removed from v.61  
changed lines
  Added in v.62
  ViewVC Help
Powered by ViewVC 1.1.21