/[lmdze]/trunk/IOIPSL/mathop.f
ViewVC logotype

Diff of /trunk/IOIPSL/mathop.f

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

trunk/libf/IOIPSL/mathop.f90 revision 32 by guez, Tue Apr 6 17:52:58 2010 UTC trunk/Sources/IOIPSL/mathop.f revision 178 by guez, Fri Mar 11 18:47:26 2016 UTC
# Line 5  MODULE mathop_m Line 5  MODULE mathop_m
5    implicit none    implicit none
6    
7    PRIVATE    PRIVATE
8    PUBLIC :: mathop    PUBLIC mathop
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 49  CONTAINS Line 56  CONTAINS
56    
57      INTEGER :: ierr      INTEGER :: ierr
58    
59      INTRINSIC SIN, COS, TAN, ASIN, ACOS, ATAN, EXP, ALOG, SQRT, ABS      !--------------------------------------------------------------------
     !---------------------------------------------------------------------  
60      ierr = 0      ierr = 0
61    
62      IF (scal >= miss_val-1.) THEN      IF (scal >= miss_val-1.) THEN
# Line 90  CONTAINS Line 96  CONTAINS
96               ierr = ma_ident_r11(nb, work_in, nb_max, work_out)               ierr = ma_ident_r11(nb, work_in, nb_max, work_out)
97            CASE DEFAULT            CASE DEFAULT
98               CALL histerr(3, "mathop", &               CALL histerr(3, "mathop", &
99                    &        'scalar variable undefined and no indexing', &          'scalar variable undefined and no indexing', &
100                    &        'but still unknown function', fun)          'but still unknown function', fun)
101            END SELECT            END SELECT
102            IF (ierr > 0) THEN            IF (ierr > 0) THEN
103               CALL histerr(3, "mathop", &               CALL histerr(3, "mathop", &
104                    &        'Error while executing a simple function', fun, ' ')          'Error while executing a simple function', fun, ' ')
105            ENDIF            ENDIF
106         ELSE         ELSE
107            SELECT CASE (fun)            SELECT CASE (fun)
108            CASE('gather')            CASE('gather')
109               ierr = ma_fugath_r11(nb, work_in, nb_index, nindex, &               ierr = ma_fugath_r11(nb, work_in, nb_index, nindex, &
110                    &                            miss_val, nb_max, work_out)                              miss_val, nb_max, work_out)
111            CASE('scatter')            CASE('scatter')
112               IF (nb_index > nb) THEN               IF (nb_index > nb) THEN
113                  work_out(1:nb_max) = miss_val                  work_out(1:nb_max) = miss_val
114                  ierr=1                  ierr=1
115               ELSE               ELSE
116                  ierr = ma_fuscat_r11(nb, work_in, nb_index, nindex, &                  ierr = ma_fuscat_r11(nb, work_in, nb_index, nindex, &
117                       &                              miss_val, nb_max, work_out)                                miss_val, nb_max, work_out)
118               ENDIF               ENDIF
119            CASE('coll')            CASE('coll')
120               ierr = ma_fucoll_r11(nb, work_in, nb_index, nindex, &               ierr = ma_fucoll_r11(nb, work_in, nb_index, nindex, nb_max, &
121                    &                            miss_val, nb_max, work_out)                    work_out)
122            CASE('fill')            CASE('fill')
123               ierr = ma_fufill_r11(nb, work_in, nb_index, nindex, &               ierr = ma_fufill_r11(nb, work_in, nb_index, nindex, nb_max, &
124                    &                            miss_val, nb_max, work_out)                    work_out)
125            CASE('undef')            CASE('undef')
126               ierr = ma_fuundef_r11(nb, work_in, nb_index, nindex, &               ierr = ma_fuundef_r11(nb, work_in, nb_index, nindex, &
127                    &                            miss_val, nb_max, work_out)                              miss_val, nb_max, work_out)
128            CASE('only')            CASE('only')
129               ierr = ma_fuonly_r11(nb, work_in, nb_index, nindex, &               ierr = ma_fuonly_r11(nb, work_in, nb_index, nindex, &
130                    &                            miss_val, nb_max, work_out)                              miss_val, nb_max, work_out)
131            CASE DEFAULT            CASE DEFAULT
132               CALL histerr(3, "mathop", &               CALL histerr(3, "mathop", &
133                    &        'scalar variable undefined and indexing', &          'scalar variable undefined and indexing', &
134                    &        'was requested but with unknown function', fun)          'was requested but with unknown function', fun)
135            END SELECT            END SELECT
136            IF (ierr > 0) THEN            IF (ierr > 0) THEN
137               CALL histerr(3, "mathop_r11", &               CALL histerr(3, "mathop_r11", &
138                    &        'Error while executing an indexing function', fun, ' ')          'Error while executing an indexing function', fun, ' ')
139            ENDIF            ENDIF
140         ENDIF         ENDIF
141      ELSE      ELSE
# Line 154  CONTAINS Line 160  CONTAINS
160            ierr = ma_power_r11(nb, work_in, scal, nb_max, work_out)            ierr = ma_power_r11(nb, work_in, scal, nb_max, work_out)
161         CASE DEFAULT         CASE DEFAULT
162            CALL histerr(3, "mathop", &            CALL histerr(3, "mathop", &
163                 &      'Unknown operation with a scalar', fun, ' ')        'Unknown operation with a scalar', fun, ' ')
164         END SELECT         END SELECT
165         IF (ierr > 0) THEN         IF (ierr > 0) THEN
166            CALL histerr(3, "mathop", &            CALL histerr(3, "mathop", &
167                 &      'Error while executing a scalar function', fun, ' ')        'Error while executing a scalar function', fun, ' ')
168         ENDIF         ENDIF
169      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  
170      !-----------------------      !-----------------------
171    END FUNCTION ma_divi_r11    END SUBROUTINE mathop_r11
172      !
173    !************************************************    SUBROUTINE mathop_r21(fun, nb, work_in, miss_val, nb_index, nindex, scal, &
174           nb_max, work_out)
175    INTEGER FUNCTION ma_power_r11(nb, x, s, nbo, y)  
176      INTEGER :: nb, nbo      ! This subroutines gives an interface to the various operations
177      REAL :: x(nb), s, y(nbo)      ! which are allowed. The interface is general enough to allow its use
178        ! for other cases.
179      INTEGER :: i  
180      !---------------------------------------------------------------------      ! INPUT
181      DO i=1, nb  
182         y(i) = x(i)**s      ! fun      : function to be applied to the vector of data
183      ENDDO      ! nb       : Length of input vector
184        ! work_in  : Input vector of data (REAL)
185      nbo = nb      ! miss_val : The value of the missing data flag (it has to be a
186      ma_power_r11 = 0      !            maximum value, in f90 : huge( a real ))
187      !-----------------------      ! nb_index : Length of index vector
188    END FUNCTION ma_power_r11      ! nindex   : Vector of indices
189        ! scal     : A scalar value for vector/scalar operations
190    !************************************************      ! 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  
191    
192      !- OUTPUT      ! OUTPUT
193    
194      !- nb_max   : Actual length of output variable      ! nb_max   : Actual length of output variable
195      !- work_out : Output vector after the operation was applied      ! work_out : Output vector after the operation was applied
196      USE errioipsl, ONLY : histerr      USE errioipsl, ONLY : histerr
197        USE mathop2, ONLY : ma_abs_r21, ma_acos_r21, ma_add_r21, ma_alog_r21, &
198             ma_asin_r21, ma_atan_r21, ma_cels_r21, ma_chs_r21, ma_cos_r21, &
199             ma_deg_r21, ma_divi_r21, ma_div_r21, ma_exp_r21, ma_fucoll_r21, &
200             ma_fufill_r21, ma_fugath_r21, ma_fumax_r21, ma_fumin_r21, &
201             ma_fuonly_r21, ma_fuscat_r21, ma_fuundef_r21, ma_ident_r21, &
202             ma_kelv_r21, ma_mult_r21, ma_power_r21, ma_rad_r21, ma_sin_r21, &
203             ma_sqrt_r21, ma_subi_r21, ma_sub_r21, ma_tan_r21
204    
205      CHARACTER(LEN=7) :: fun      CHARACTER(LEN=7) :: fun
206      INTEGER :: nb(2), nb_max, nb_index      INTEGER :: nb(2), nb_max, nb_index
# Line 789  CONTAINS Line 210  CONTAINS
210    
211      INTEGER :: ierr      INTEGER :: ierr
212    
213      INTRINSIC SIN, COS, TAN, ASIN, ACOS, ATAN, EXP, ALOG, SQRT, ABS      !--------------------------------------------------------------------
     !---------------------------------------------------------------------  
214      ierr = 0      ierr = 0
215    
216      IF (scal >= miss_val-1.) THEN      IF (scal >= miss_val-1.) THEN
# Line 830  CONTAINS Line 250  CONTAINS
250               ierr = ma_ident_r21(nb, work_in, nb_max, work_out)               ierr = ma_ident_r21(nb, work_in, nb_max, work_out)
251            CASE DEFAULT            CASE DEFAULT
252               CALL histerr(3, "mathop", &               CALL histerr(3, "mathop", &
253                    &        'scalar variable undefined and no indexing', &          'scalar variable undefined and no indexing', &
254                    &        'but still unknown function', fun)          'but still unknown function', fun)
255            END SELECT            END SELECT
256            IF (ierr > 0) THEN            IF (ierr > 0) THEN
257               CALL histerr(3, "mathop", &               CALL histerr(3, "mathop", &
258                    &        'Error while executing a simple function', fun, ' ')          'Error while executing a simple function', fun, ' ')
259            ENDIF            ENDIF
260         ELSE         ELSE
261            SELECT CASE (fun)            SELECT CASE (fun)
262            CASE('gather')            CASE('gather')
263               ierr = ma_fugath_r21(nb, work_in, nb_index, nindex, &               ierr = ma_fugath_r21(nb, work_in, nb_index, nindex, &
264                    &                            miss_val, nb_max, work_out)                              miss_val, nb_max, work_out)
265            CASE('scatter')            CASE('scatter')
266               IF (nb_index > (nb(1)*nb(2)) ) THEN               IF (nb_index > (nb(1)*nb(2)) ) THEN
267                  work_out(1:nb_max) = miss_val                  work_out(1:nb_max) = miss_val
268                  ierr=1                  ierr=1
269               ELSE               ELSE
270                  ierr = ma_fuscat_r21(nb, work_in, nb_index, nindex, &                  ierr = ma_fuscat_r21(nb, work_in, nb_index, nindex, &
271                       &                             miss_val, nb_max, work_out)                               miss_val, nb_max, work_out)
272               ENDIF               ENDIF
273            CASE('coll')            CASE('coll')
274               ierr = ma_fucoll_r21(nb, work_in, nb_index, nindex, &               ierr = ma_fucoll_r21(nb, work_in, nb_index, nindex, nb_max, &
275                    &                           miss_val, nb_max, work_out)                    work_out)
276            CASE('fill')            CASE('fill')
277               ierr = ma_fufill_r21(nb, work_in, nb_index, nindex, &               ierr = ma_fufill_r21(nb, work_in, nb_index, nindex, nb_max, &
278                    &                           miss_val, nb_max, work_out)                    work_out)
279            CASE('undef')            CASE('undef')
280               ierr = ma_fuundef_r21(nb, work_in, nb_index, nindex, &               ierr = ma_fuundef_r21(nb, work_in, nb_index, nindex, &
281                    &                           miss_val, nb_max, work_out)                             miss_val, nb_max, work_out)
282            CASE('only')            CASE('only')
283               ierr = ma_fuonly_r21(nb, work_in, nb_index, nindex, &               ierr = ma_fuonly_r21(nb, work_in, nb_index, nindex, &
284                    &                           miss_val, nb_max, work_out)                             miss_val, nb_max, work_out)
285            CASE DEFAULT            CASE DEFAULT
286               CALL histerr(3, "mathop", &               CALL histerr(3, "mathop", &
287                    &        'scalar variable undefined and indexing', &          'scalar variable undefined and indexing', &
288                    &        'was requested but with unknown function', fun)          'was requested but with unknown function', fun)
289            END SELECT            END SELECT
290            IF (ierr > 0) THEN            IF (ierr > 0) THEN
291               CALL histerr(3, "mathop_r21", &               CALL histerr(3, "mathop_r21", &
292                    &        'Error while executing an indexing function', fun, ' ')          'Error while executing an indexing function', fun, ' ')
293            ENDIF            ENDIF
294         ENDIF         ENDIF
295      ELSE      ELSE
# Line 894  CONTAINS Line 314  CONTAINS
314            ierr = ma_power_r21(nb, work_in, scal, nb_max, work_out)            ierr = ma_power_r21(nb, work_in, scal, nb_max, work_out)
315         CASE DEFAULT         CASE DEFAULT
316            CALL histerr(3, "mathop", &            CALL histerr(3, "mathop", &
317                 &      'Unknown operation with a scalar', fun, ' ')        'Unknown operation with a scalar', fun, ' ')
318         END SELECT         END SELECT
319         IF (ierr > 0) THEN         IF (ierr > 0) THEN
320            CALL histerr(3, "mathop", &            CALL histerr(3, "mathop", &
321                 &      'Error while executing a scalar function', fun, ' ')        'Error while executing a scalar function', fun, ' ')
322         ENDIF         ENDIF
323      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  
324      !-----------------------      !-----------------------
325    END FUNCTION ma_kelv_r21    END SUBROUTINE mathop_r21
326      !
327    !************************************************    SUBROUTINE mathop_r31(fun, nb, work_in, miss_val, nb_index, nindex, scal, &
328           nb_max, work_out)
329    INTEGER FUNCTION ma_deg_r21(nb, x, nbo, y)  
330      INTEGER :: nb(2), nbo, i, j, ij      ! This subroutines gives an interface to the various operations
331      REAL :: x(nb(1), nb(2)), y(nbo)      ! which are allowed. The interface is general enough to allow its use
332      !---------------------------------------------------------------------      ! for other cases.
333      ij = 0  
334      DO j=1, nb(2)      ! INPUT
335         DO i=1, nb(1)  
336            ij = ij+1      ! fun      : function to be applied to the vector of data
337            y(ij) = x(i, j)*57.29577951      ! nb       : Length of input vector
338         ENDDO      ! work_in  : Input vector of data (REAL)
339      ENDDO      ! miss_val : The value of the missing data flag (it has to be a
340        !            maximum value, in f90 : huge( a real ))
341      nbo = nb(1)*nb(2)      ! nb_index : Length of index vector
342      ma_deg_r21 = 0      ! nindex   : Vector of indices
343      !----------------------      ! scal     : A scalar value for vector/scalar operations
344    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  
345    
346      !- OUTPUT      ! OUTPUT
347    
348      !- nb_max   : Actual length of output variable      ! nb_max   : Actual length of output variable
349      !- work_out : Output vector after the operation was applied      ! work_out : Output vector after the operation was applied
350      USE errioipsl, ONLY : histerr      USE errioipsl, ONLY : histerr
351        USE mathop2, ONLY : ma_abs_r31, ma_acos_r31, ma_add_r31, ma_alog_r31, &
352             ma_asin_r31, ma_atan_r31, ma_cels_r31, ma_chs_r31, ma_cos_r31, &
353             ma_deg_r31, ma_divi_r31, ma_div_r31, ma_exp_r31, ma_fucoll_r31, &
354             ma_fufill_r31, ma_fugath_r31, ma_fumax_r31, ma_fumin_r31, &
355             ma_fuonly_r31, ma_fuscat_r31, ma_fuundef_r31, ma_ident_r31, &
356             ma_kelv_r31, ma_mult_r31, ma_power_r31, ma_rad_r31, ma_sin_r31, &
357             ma_sqrt_r31, ma_subi_r31, ma_sub_r31, ma_tan_r31
358    
359      CHARACTER(LEN=7) :: fun      CHARACTER(LEN=7) :: fun
360      INTEGER :: nb(3), nb_max, nb_index      INTEGER :: nb(3), nb_max, nb_index
# Line 1639  CONTAINS Line 364  CONTAINS
364    
365      INTEGER :: ierr      INTEGER :: ierr
366    
367      INTRINSIC SIN, COS, TAN, ASIN, ACOS, ATAN, EXP, ALOG, SQRT, ABS      !--------------------------------------------------------------------
     !---------------------------------------------------------------------  
368      ierr = 0      ierr = 0
369    
370      IF (scal >= miss_val-1.) THEN      IF (scal >= miss_val-1.) THEN
# Line 1680  CONTAINS Line 404  CONTAINS
404               ierr = ma_ident_r31(nb, work_in, nb_max, work_out)               ierr = ma_ident_r31(nb, work_in, nb_max, work_out)
405            CASE DEFAULT            CASE DEFAULT
406               CALL histerr(3, "mathop", &               CALL histerr(3, "mathop", &
407                    &        'scalar variable undefined and no indexing', &                    'scalar variable undefined and no indexing', &
408                    &        'but still unknown function', fun)                    'but still unknown function', fun)
409            END SELECT            END SELECT
410            IF (ierr > 0) THEN            IF (ierr > 0) THEN
411               CALL histerr(3, "mathop", &               CALL histerr(3, "mathop", &
412                    &        'Error while executing a simple function', fun, ' ')                    'Error while executing a simple function', fun, ' ')
413            ENDIF            ENDIF
414         ELSE         ELSE
415            SELECT CASE (fun)            SELECT CASE (fun)
416            CASE('gather')            CASE('gather')
417               ierr = ma_fugath_r31(nb, work_in, nb_index, nindex, &               ierr = ma_fugath_r31(nb, work_in, nb_index, nindex, &
418                    &                           miss_val, nb_max, work_out)                    miss_val, nb_max, work_out)
419            CASE('scatter')            CASE('scatter')
420               IF (nb_index > (nb(1)*nb(2)*nb(3))) THEN               IF (nb_index > (nb(1)*nb(2)*nb(3))) THEN
421                  work_out(1:nb_max) = miss_val                  work_out(1:nb_max) = miss_val
422                  ierr=1                  ierr=1
423               ELSE               ELSE
424                  ierr = ma_fuscat_r31(nb, work_in, nb_index, nindex, &                  ierr = ma_fuscat_r31(nb, work_in, nb_index, nindex, &
425                       &                             miss_val, nb_max, work_out)                       miss_val, nb_max, work_out)
426               ENDIF               ENDIF
427            CASE('coll')            CASE('coll')
428               ierr = ma_fucoll_r31(nb, work_in, nb_index, nindex, &               ierr = ma_fucoll_r31(nb, work_in, nb_index, nindex, nb_max, &
429                    &                           miss_val, nb_max, work_out)                    work_out)
430            CASE('fill')            CASE('fill')
431               ierr = ma_fufill_r31(nb, work_in, nb_index, nindex, &               ierr = ma_fufill_r31(nb, work_in, nb_index, nindex, nb_max, &
432                    &                           miss_val, nb_max, work_out)                    work_out)
433            CASE('undef')            CASE('undef')
434               ierr = ma_fuundef_r31(nb, work_in, nb_index, nindex, &               ierr = ma_fuundef_r31(nb, work_in, nb_index, nindex, &
435                    &                           miss_val, nb_max, work_out)                    miss_val, nb_max, work_out)
436            CASE('only')            CASE('only')
437               ierr = ma_fuonly_r31(nb, work_in, nb_index, nindex, &               ierr = ma_fuonly_r31(nb, work_in, nb_index, nindex, &
438                    &                           miss_val, nb_max, work_out)                    miss_val, nb_max, work_out)
439            CASE DEFAULT            CASE DEFAULT
440               CALL histerr(3, "mathop", &               CALL histerr(3, "mathop", &
441                    &        'scalar variable undefined and indexing', &                    'scalar variable undefined and indexing', &
442                    &        'was requested but with unknown function', fun)                    'was requested but with unknown function', fun)
443            END SELECT            END SELECT
444            IF (ierr > 0) THEN            IF (ierr > 0) THEN
445               CALL histerr(3, "mathop_r31", &               CALL histerr(3, "mathop_r31", &
446                    &        'Error while executing an indexing function', fun, ' ')                    'Error while executing an indexing function', fun, ' ')
447            ENDIF            ENDIF
448         ENDIF         ENDIF
449      ELSE      ELSE
# Line 1744  CONTAINS Line 468  CONTAINS
468            ierr = ma_power_r31(nb, work_in, scal, nb_max, work_out)            ierr = ma_power_r31(nb, work_in, scal, nb_max, work_out)
469         CASE DEFAULT         CASE DEFAULT
470            CALL histerr(3, "mathop", &            CALL histerr(3, "mathop", &
471                 &      'Unknown operation with a scalar', fun, ' ')                 'Unknown operation with a scalar', fun, ' ')
472         END SELECT         END SELECT
473         IF (ierr > 0) THEN         IF (ierr > 0) THEN
474            CALL histerr(3, "mathop", &            CALL histerr(3, "mathop", &
475                 &      'Error while executing a scalar function', fun, ' ')                 'Error while executing a scalar function', fun, ' ')
476         ENDIF         ENDIF
477      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  
   
   !************************************************  
478    
479    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  
480    
481  END MODULE mathop_m  END MODULE mathop_m
Legend:
Removed from v.32  
changed lines
  Added in v.178
  ViewVC Help
Powered by ViewVC 1.1.21