ns5001k-admin/Nemo-NAA10km_nemo_4.2.2
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Nemo-NAA10km_nemo_4.2.2/tools/REBUILD/BLD/ppsrc/nemo/mathelp.f90
Robinson Hordoir 4fa5daa6c3 First commit of Nemo-NAA10km_nemo_4.2.2
2025-09-05 09:42:21 +02:00

3134 lines
80 KiB
Fortran
Raw Permalink Blame History

MODULE mathelp
!-
!$Id: mathelp.f90 2281 2010-10-15 14:21:13Z smasson $
!-
! This software is governed by the CeCILL license
! See IOIPSL/IOIPSL_License_CeCILL.txt
!---------------------------------------------------------------------
USE errioipsl,ONLY : ipslerr
USE stringop
!-
PRIVATE
PUBLIC :: mathop,moycum,buildop
!-
INTERFACE mathop
MODULE PROCEDURE mathop_r11,mathop_r21,mathop_r31
END INTERFACE
!-
!- Variables used to detect and identify the operations
!-
CHARACTER(LEN=80),SAVE :: &
& seps='( ) , + - / * ^', ops = '+ - * / ^', mima = 'min max'
CHARACTER(LEN=250),SAVE :: &
& funcs = 'sin cos tan asin acos atan exp log sqrt chs abs '&
& //'cels kelv deg rad gather scatter fill coll undef only ident'
CHARACTER(LEN=120),SAVE :: &
& indexfu = 'gather, scatter, fill, coll, undef, only'
!---------------------------------------------------------------------
CONTAINS
!===
SUBROUTINE buildop (c_str,ex_topps,topp,fill_val,opps,scal,nbops)
!---------------------------------------------------------------------
!- This subroutine decomposes the input string in the elementary
!- functions which need to be applied to the vector of data.
!- This vector is represented by X in the string.
!- This subroutine is the driver of the decomposition and gets
!- the time operation but then call decoop for the other operations
!-
!- INPUT
!-
!- c_str : String containing the operations
!- ex_toops : Time operations that can be expected within the string
!- fill_val :
!-
!- OUTPUT
!-
!- topp : Time operation
!- opps :
!- scal :
!- nbops :
!---------------------------------------------------------------------
IMPLICIT NONE
!-
CHARACTER(LEN=*),INTENT(IN) :: c_str,ex_topps
CHARACTER(LEN=*),INTENT(OUT) :: topp
CHARACTER(LEN=*),DIMENSION(:),INTENT(OUT) :: opps
REAL,INTENT(IN) :: fill_val
REAL,DIMENSION(:),INTENT(OUT) :: scal
INTEGER,INTENT(OUT) :: nbops
!-
CHARACTER(LEN=LEN(c_str)) :: str,new_str
INTEGER :: leng,ind_opb,ind_clb
!-
LOGICAL :: check = .FALSE.
!---------------------------------------------------------------------
IF (check) WRITE(*,*) 'buildop : Some preliminary cleaning'
!-
str = c_str
leng = LEN_TRIM(str)
IF ( str(1:1) == '(' .AND. str(leng:leng) == ')' ) THEN
str = str(2:leng-1)
leng = leng-2
ENDIF
!-
IF (check) &
& WRITE(*,*) 'buildop : Starting to test the various options'
!-
IF (leng <= 5 .AND. INDEX(ex_topps,str(1:leng)) > 0) THEN
IF (check) WRITE(*,*) 'buildop : Time operation only'
nbops = 0
topp = str(1:leng)
ELSE
IF (check) THEN
WRITE(*,*) 'buildop : Time operation and something else'
ENDIF
!--
ind_opb = INDEX(str(1:leng),'(')
IF (ind_opb > 0) THEN
IF (INDEX(ex_topps,str(1:ind_opb-1)) > 0) THEN
IF (check) THEN
WRITE(*,'(2a)') &
& ' buildop : Extract time operation from : ',str
ENDIF
topp = str(1:ind_opb-1)
ind_clb = INDEX(str(1:leng),')',BACK=.TRUE.)
new_str = str(ind_opb+1:ind_clb-1)
IF (check) THEN
WRITE(*,'(2a,2I3)') &
& ' buildop : Call decoop ',new_str,ind_opb,ind_clb
ENDIF
CALL decoop (new_str,fill_val,opps,scal,nbops)
ELSE
CALL ipslerr(3,'buildop', &
& 'time operation does not exist',str(1:ind_opb-1),' ')
ENDIF
ELSE
CALL ipslerr(3,'buildop', &
& 'some long operation exists but wihout parenthesis', &
& str(1:leng),' ')
ENDIF
ENDIF
!-
IF (check) THEN
DO leng=1,nbops
WRITE(*,*) &
& 'buildop : i -- opps, scal : ',leng,opps(leng),scal(leng)
ENDDO
ENDIF
!---------------------
END SUBROUTINE buildop
!===
SUBROUTINE decoop (pstr,fill_val,opps,scal,nbops)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
CHARACTER(LEN=*),INTENT(IN) :: pstr
REAL,INTENT(IN) :: fill_val
CHARACTER(LEN=*),DIMENSION(:),INTENT(OUT) :: opps
REAL,DIMENSION(:),INTENT(OUT) :: scal
INTEGER,INTENT(OUT) :: nbops
!-
CHARACTER(LEN=1),DIMENSION(2) :: f_char,s_char
INTEGER,DIMENSION(2) :: f_pos,s_pos
CHARACTER(LEN=20) :: opp_str,scal_str
CHARACTER(LEN=LEN(pstr)) :: str
INTEGER :: nbsep,nbops_max,xpos,leng,ppos,epos,int_tmp
CHARACTER(LEN=3) :: tl,dl
CHARACTER(LEN=10) :: fmt
!-
LOGICAL :: check = .FALSE.,prio
!---------------------------------------------------------------------
IF (check) WRITE(*,'(2A)') ' decoop : Incoming string : ',pstr
!-
str = pstr; nbops = 0;
!-
CALL findsep (str,nbsep,f_char,f_pos,s_char,s_pos)
IF (check) WRITE(*,*) 'decoop : Out of findsep',nbsep
!-
nbops_max = min(SIZE(opps),SIZE(scal))
!-
DO WHILE (nbsep > 0)
IF (nbops >= nbops_max) THEN
CALL ipslerr(3,'decoop','Expression too complex',TRIM(str),' ')
ENDIF
!--
xpos = INDEX(str,'X')
leng = LEN_TRIM(str)
nbops = nbops+1
!--
IF (check) THEN
WRITE(*,*) 'decoop : str -> ',TRIM(str)
WRITE(*,*) 'decoop : nbops -> ',nbops
WRITE(*,*) s_char(1),'-',f_char(1),'|',f_char(2),'-',s_char(2)
WRITE(*,*) s_pos(1),'-',f_pos(1),'|',f_pos(2),'-',s_pos(2)
ENDIF
!---
!-- Start the analysis of the syntax. 3 types of constructs
!-- are recognized. They are scanned sequentialy
!---
IF (nbsep == 1) THEN
IF (check) WRITE(*,*) 'decoop : Only one operation'
IF (INDEX(ops,f_char(1)) > 0) THEN
!------ Type : scal+X
IF (f_char(1) == '-' .OR. f_char(1) == '/') THEN
opp_str = f_char(1)//'I'
ELSE
opp_str = f_char(1)
ENDIF
scal_str = str(s_pos(1)+1:f_pos(1)-1)
str = 'X'
ELSE IF (INDEX(ops,f_char(2)) > 0) THEN
!------ Type : X+scal
opp_str = f_char(2)
scal_str = str(f_pos(2)+1:s_pos(2)-1)
str = 'X'
ELSE
CALL ipslerr(3,'decoop', &
& 'Unknown operations of type X+scal',f_char(1),pstr)
ENDIF
ELSE
IF (check) WRITE(*,*) 'decoop : More complex operation'
IF ( f_char(1) == '(' .AND. f_char(2) == ')' ) THEN
!------ Type : sin(X)
opp_str = str(s_pos(1)+1:f_pos(1)-1)
scal_str = '?'
str = str(1:s_pos(1))//'X'//str(f_pos(2)+1:leng)
ELSE IF ( (f_char(1) == '(' .AND. f_char(2) == ',')&
& .OR.(f_char(1) == ',' .AND. f_char(2) == ')')) THEN
!------ Type : max(X,scal) or max(scal,X)
IF (f_char(1) == '(' .AND. s_char(2) == ')') THEN
!-------- Type : max(X,scal)
opp_str = str(f_pos(1)-3:f_pos(1)-1)
scal_str = str(f_pos(2)+1:s_pos(2)-1)
str = str(1:f_pos(1)-4)//'X'//str(s_pos(2)+1:leng)
ELSE IF (f_char(1) == ',' .AND. s_char(1) == '(') THEN
!-------- Type : max(scal,X)
opp_str = str(s_pos(1)-3:s_pos(1)-1)
scal_str = str(s_pos(1)+1:f_pos(1)-1)
str = str(1:s_pos(1)-4)//'X'//str(f_pos(2)+1:leng)
ELSE
CALL ipslerr(3,'decoop','Syntax error 1',str,' ')
ENDIF
ELSE
prio = (f_char(2) == '*').OR.(f_char(2) == '^')
IF ( (INDEX(ops,f_char(1)) > 0) &
& .AND.(xpos-f_pos(1) == 1).AND.(.NOT.prio) ) THEN
!-------- Type : ... scal+X ...
IF (f_char(1) == '-' .OR. f_char(1) == '/') THEN
opp_str = f_char(1)//'I'
ELSE
opp_str = f_char(1)
ENDIF
scal_str = str(s_pos(1)+1:f_pos(1)-1)
str = str(1:s_pos(1))//'X'//str(f_pos(1)+2:leng)
ELSE IF ( (INDEX(ops,f_char(2)) > 0) &
& .AND.(f_pos(2)-xpos == 1) ) THEN
!-------- Type : ... X+scal ...
opp_str = f_char(2)
scal_str = str(f_pos(2)+1:s_pos(2)-1)
str = str(1:f_pos(2)-2)//'X'//str(s_pos(2):leng)
ELSE
CALL ipslerr(3,'decoop','Syntax error 2',str,' ')
ENDIF
ENDIF
ENDIF
!---
IF (check) WRITE(*,*) 'decoop : Finished syntax,str = ',TRIM(str)
!---
!-- Now that the different components of the operation are identified
!-- we transform them into what is going to be used in the program
!---
IF (INDEX(scal_str,'?') > 0) THEN
IF (INDEX(funcs,opp_str(1:LEN_TRIM(opp_str))) > 0) THEN
opps(nbops) = opp_str(1:LEN_TRIM(opp_str))
scal(nbops) = fill_val
ELSE
CALL ipslerr(3,'decoop', &
& 'Unknown function',opp_str(1:LEN_TRIM(opp_str)),' ')
ENDIF
ELSE
leng = LEN_TRIM(opp_str)
IF (INDEX(mima,opp_str(1:leng)) > 0) THEN
opps(nbops) = 'fu'//opp_str(1:leng)
ELSE
IF (INDEX(opp_str(1:leng),'+') > 0) THEN
opps(nbops) = 'add'
ELSE IF (INDEX(opp_str(1:leng),'-I') > 0) THEN
opps(nbops) = 'subi'
ELSE IF (INDEX(opp_str(1:leng),'-') > 0) THEN
opps(nbops) = 'sub'
ELSE IF (INDEX(opp_str(1:leng),'*') > 0) THEN
opps(nbops) = 'mult'
ELSE IF (INDEX(opp_str(1:leng),'/') > 0) THEN
opps(nbops) = 'div'
ELSE IF (INDEX(opp_str(1:leng),'/I') > 0) THEN
opps(nbops) = 'divi'
ELSE IF (INDEX(opp_str(1:leng),'^') > 0) THEN
opps(nbops) = 'power'
ELSE
CALL ipslerr(3,'decoop', &
& 'Unknown operation',opp_str(1:leng),' ')
ENDIF
ENDIF
!-----
leng = LEN_TRIM(scal_str)
ppos = INDEX(scal_str,'.')
epos = INDEX(scal_str,'e')
IF (epos == 0) epos = INDEX(scal_str,'E')
!-----
!---- Try to catch a few errors
!-----
IF (INDEX(ops,scal_str) > 0) THEN
CALL ipslerr(3,'decoop', &
& 'Strange scalar you have here ',scal_str,pstr)
ENDIF
IF (epos > 0) THEN
WRITE(tl,'(I3.3)') leng
WRITE(dl,'(I3.3)') epos-ppos-1
fmt='(e'//tl//'.'//dl//')'
READ(scal_str,fmt) scal(nbops)
ELSE IF (ppos > 0) THEN
WRITE(tl,'(I3.3)') leng
WRITE(dl,'(I3.3)') leng-ppos
fmt='(f'//tl//'.'//dl//')'
READ(scal_str,fmt) scal(nbops)
ELSE
WRITE(tl,'(I3.3)') leng
fmt = '(I'//tl//')'
READ(scal_str,fmt) int_tmp
scal(nbops) = REAL(int_tmp)
ENDIF
ENDIF
IF (check) WRITE(*,*) 'decoop : Finished interpretation'
CALL findsep(str,nbsep,f_char,f_pos,s_char,s_pos)
ENDDO
!--------------------
END SUBROUTINE decoop
!===
SUBROUTINE findsep (str,nbsep,f_char,f_pos,s_char,s_pos)
!---------------------------------------------------------------------
!- Subroutine finds all separators in a given string
!- It returns the following information about str :
!- f_char : The first separation character
!- (1 for before and 2 for after)
!- f_pos : The position of the first separator
!- s_char : The second separation character
!- (1 for before and 2 for after)
!- s_pos : The position of the second separator
!---------------------------------------------------------------------
IMPLICIT NONE
!-
CHARACTER(LEN=*),INTENT(INOUT) :: str
INTEGER :: nbsep
CHARACTER(LEN=1),DIMENSION(2) :: f_char,s_char
INTEGER,DIMENSION(2) :: f_pos,s_pos
!-
CHARACTER(LEN=10) :: str_tmp
LOGICAL :: f_found,s_found
INTEGER :: ind,xpos,leng,i
!-
LOGICAL :: check = .FALSE.
!---------------------------------------------------------------------
IF (check) WRITE(*,*) 'findsep : call cleanstr: ',TRIM(str)
!-
CALL cleanstr(str)
!-
IF (check) WRITE(*,*) 'findsep : out of cleanstr: ',TRIM(str)
!-
xpos = INDEX(str,'X')
leng = LEN_TRIM(str)
!-
f_pos(1:2) = (/ 0,leng+1 /)
f_char(1:2) = (/ '?','?' /)
s_pos(1:2) = (/ 0,leng+1 /)
s_char(1:2) = (/ '?','?' /)
!-
nbsep = 0
!-
f_found = .FALSE.
s_found = .FALSE.
IF (xpos > 1) THEN
DO i=xpos-1,1,-1
ind = INDEX(seps,str(i:i))
IF (ind > 0) THEN
IF (.NOT.f_found) THEN
f_char(1) = str(i:i)
f_pos(1) = i
nbsep = nbsep+1
f_found = .TRUE.
ELSE IF (.NOT.s_found) THEN
s_char(1) = str(i:i)
s_pos(1) = i
nbsep = nbsep+1
s_found = .TRUE.
ENDIF
ENDIF
ENDDO
ENDIF
!-
f_found = .FALSE.
s_found = .FALSE.
IF (xpos < leng) THEN
DO i=xpos+1,leng
ind = INDEX(seps,str(i:i))
IF (ind > 0) THEN
IF (.NOT.f_found) THEN
f_char(2) = str(i:i)
f_pos(2) = i
nbsep = nbsep+1
f_found = .TRUE.
ELSE IF (.NOT.s_found) THEN
s_char(2) = str(i:i)
s_pos(2) = i
nbsep = nbsep+1
s_found = .TRUE.
ENDIF
ENDIF
ENDDO
ENDIF
!-
IF (nbsep > 4) THEN
WRITE(str_tmp,'("number :",I3)') nbsep
CALL ipslerr(3,'findsep', &
& 'How can I find that many separators',str_tmp,TRIM(str))
ENDIF
!-
IF (check) WRITE(*,*) 'Finished findsep : ',nbsep,leng
!---------------------
END SUBROUTINE findsep
!===
SUBROUTINE cleanstr(str)
!---------------------------------------------------------------------
!- We clean up the string by taking out the extra () and puting
!- everything in lower case except for the X describing the variable
!---------------------------------------------------------------------
IMPLICIT NONE
!-
CHARACTER(LEN=*),INTENT(INOUT) :: str
!-
INTEGER :: ind,leng,ic,it
LOGICAL :: check = .FALSE.
!---------------------------------------------------------------------
leng = LEN_TRIM(str)
CALL strlowercase(str)
!-
ind = INDEX(str,'x')
IF (check) THEN
WRITE (*,*) 'cleanstr 1.0 : ind = ',ind, &
& ' str = ',str(1:leng),'---'
ENDIF
!-
! If the character before the x is not a letter then we can assume
! that it is the variable and promote it to a capital letter
!-
DO WHILE (ind > 0)
ic = 0
IF (ind > 1) ic = IACHAR(str(ind-1:ind-1))
IF (ic < 97 .OR. ic > 122) THEN
str(ind:ind) = 'X'
ENDIF
it = INDEX(str(ind+1:leng),'x')
IF (it > 0) THEN
ind = ind+it
ELSE
ind = it
ENDIF
ENDDO
!-
IF (check) WRITE (*,*) 'cleanstr 2.0 : str = ',str(1:leng),'---'
!-
IF ( str(1:1) == '(' .AND. str(leng:leng) == ')' ) THEN
str = str(2:leng-1)
ENDIF
!-
IF (check) WRITE (*,*) 'cleanstr 3.0 : str = ',str(1:leng),'---'
!-
leng = LEN_TRIM(str)
ind = INDEX(str,'((X))')
IF (ind > 0) THEN
str=str(1:ind-1)//'(X)'//str(ind+5:leng)//' '
ENDIF
!-
IF (check) WRITE (*,*) 'cleanstr 4.0 : str = ',str(1:leng),'---'
!-
leng = LEN_TRIM(str)
ind = INDEX(str,'(X)')
IF (ind > 0 .AND. ind+3 < leng) THEN
IF ( (INDEX(seps,str(ind-1:ind-1)) > 0) &
& .AND. (INDEX(seps,str(ind+3:ind+3)) > 0) ) THEN
str=str(1:ind-1)//'X'//str(ind+3:leng)//' '
ENDIF
ENDIF
!-
IF (check) WRITE (*,*) 'cleanstr 5.0 : str = ',str(1:leng),'---'
!-
leng = LEN_TRIM(str)
ind = INDEX(str(1:leng),' ')
DO WHILE (ind > 0)
str=str(1:ind-1)//str(ind+1:leng)//' '
leng = LEN_TRIM(str)
ind = INDEX(str(1:leng),' ')
ENDDO
!-
IF (check) WRITE (*,*) 'cleanstr 6.0 : str = ',str(1:leng),'---'
!----------------------
END SUBROUTINE cleanstr
!===
!===
SUBROUTINE mathop_r11 &
& (fun,nb,work_in,miss_val,nb_index,nindex,scal,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
!-
!- OUTPUT
!-
!- nb_max : Actual length of output variable
!- work_out : Output vector after the operation was applied
!---------------------------------------------------------------------
IMPLICIT NONE
!-
CHARACTER(LEN=7) :: fun
INTEGER :: nb,nb_max,nb_index
INTEGER :: nindex(nb_index)
REAL :: work_in(nb),scal,miss_val
REAL :: work_out(nb_max)
!-
INTEGER :: ierr
!-
INTRINSIC SIN,COS,TAN,ASIN,ACOS,ATAN,EXP,LOG,SQRT,ABS
!---------------------------------------------------------------------
ierr = 0
!-
IF (scal >= miss_val-1.) THEN
IF (INDEX(indexfu,fun(1:LEN_TRIM(fun))) == 0) THEN
SELECT CASE (fun)
CASE('sin')
ierr = ma_sin_r11(nb,work_in,nb_max,work_out)
CASE('cos')
ierr = ma_cos_r11(nb,work_in,nb_max,work_out)
CASE('tan')
ierr = ma_tan_r11(nb,work_in,nb_max,work_out)
CASE('asin')
ierr = ma_asin_r11(nb,work_in,nb_max,work_out)
CASE('acos')
ierr = ma_acos_r11(nb,work_in,nb_max,work_out)
CASE('atan')
ierr = ma_atan_r11(nb,work_in,nb_max,work_out)
CASE('exp')
ierr = ma_exp_r11(nb,work_in,nb_max,work_out)
CASE('log')
ierr = ma_log_r11(nb,work_in,nb_max,work_out)
CASE('sqrt')
ierr = ma_sqrt_r11(nb,work_in,nb_max,work_out)
CASE('chs')
ierr = ma_chs_r11(nb,work_in,nb_max,work_out)
CASE('abs')
ierr = ma_abs_r11(nb,work_in,nb_max,work_out)
CASE('cels')
ierr = ma_cels_r11(nb,work_in,nb_max,work_out)
CASE('kelv')
ierr = ma_kelv_r11(nb,work_in,nb_max,work_out)
CASE('deg')
ierr = ma_deg_r11(nb,work_in,nb_max,work_out)
CASE('rad')
ierr = ma_rad_r11(nb,work_in,nb_max,work_out)
CASE('ident')
ierr = ma_ident_r11(nb,work_in,nb_max,work_out)
CASE DEFAULT
CALL ipslerr(3,"mathop", &
& 'scalar variable undefined and no indexing', &
& 'but still unknown function',fun)
END SELECT
IF (ierr > 0) THEN
CALL ipslerr(3,"mathop", &
& 'Error while executing a simple function',fun,' ')
ENDIF
ELSE
SELECT CASE (fun)
CASE('gather')
ierr = ma_fugath_r11(nb,work_in,nb_index,nindex, &
& miss_val,nb_max,work_out)
CASE('scatter')
IF (nb_index > nb) THEN
work_out(1:nb_max) = miss_val
ierr=1
ELSE
ierr = ma_fuscat_r11(nb,work_in,nb_index,nindex, &
& miss_val,nb_max,work_out)
ENDIF
CASE('coll')
ierr = ma_fucoll_r11(nb,work_in,nb_index,nindex, &
& miss_val,nb_max,work_out)
CASE('fill')
ierr = ma_fufill_r11(nb,work_in,nb_index,nindex, &
& miss_val,nb_max,work_out)
CASE('undef')
ierr = ma_fuundef_r11(nb,work_in,nb_index,nindex, &
& miss_val,nb_max,work_out)
CASE('only')
ierr = ma_fuonly_r11(nb,work_in,nb_index,nindex, &
& miss_val,nb_max,work_out)
CASE DEFAULT
CALL ipslerr(3,"mathop", &
& 'scalar variable undefined and indexing',&
& 'was requested but with unknown function',fun)
END SELECT
IF (ierr > 0) THEN
CALL ipslerr(3,"mathop_r11", &
& 'Error while executing an indexing function',fun,' ')
ENDIF
ENDIF
ELSE
SELECT CASE (fun)
CASE('fumin')
ierr = ma_fumin_r11(nb,work_in,scal,nb_max,work_out)
CASE('fumax')
ierr = ma_fumax_r11(nb,work_in,scal,nb_max,work_out)
CASE('add')
ierr = ma_add_r11(nb,work_in,scal,nb_max,work_out)
CASE('subi')
ierr = ma_subi_r11(nb,work_in,scal,nb_max,work_out)
CASE('sub')
ierr = ma_sub_r11(nb,work_in,scal,nb_max,work_out)
CASE('mult')
ierr = ma_mult_r11(nb,work_in,scal,nb_max,work_out)
CASE('div')
ierr = ma_div_r11(nb,work_in,scal,nb_max,work_out)
CASE('divi')
ierr = ma_divi_r11(nb,work_in,scal,nb_max,work_out)
CASE('power')
ierr = ma_power_r11(nb,work_in,scal,nb_max,work_out)
CASE DEFAULT
CALL ipslerr(3,"mathop", &
& 'Unknown operation with a scalar',fun,' ')
END SELECT
IF (ierr > 0) THEN
CALL ipslerr(3,"mathop", &
& 'Error while executing a scalar function',fun,' ')
ENDIF
ENDIF
!------------------------
END SUBROUTINE mathop_r11
!-
!=== FUNCTIONS (only one argument)
!-
INTEGER FUNCTION ma_sin_r11(nb,x,nbo,y)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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_log_r11(nb,x,nbo,y)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
INTEGER :: nb,nbo,i
REAL :: x(nb),y(nbo)
!---------------------------------------------------------------------
DO i=1,nb
y(i) = log(x(i))
ENDDO
!-
nbo = nb
ma_log_r11 = 0
!----------------------
END FUNCTION ma_log_r11
!===
INTEGER FUNCTION ma_sqrt_r11(nb,x,nbo,y)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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
!-----------------------
END FUNCTION ma_divi_r11
!===
INTEGER FUNCTION ma_power_r11(nb,x,s,nbo,y)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
INTEGER :: nb,nbo
REAL :: x(nb),s,y(nbo)
!-
INTEGER :: i
!---------------------------------------------------------------------
DO i=1,nb
y(i) = x(i)**s
ENDDO
!-
nbo = nb
ma_power_r11 = 0
!-----------------------
END FUNCTION ma_power_r11
!===
INTEGER FUNCTION ma_fumin_r11(nb,x,s,nbo,y)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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
!-
!- OUTPUT
!-
!- nb_max : Actual length of output variable
!- work_out : Output vector after the operation was applied
!---------------------------------------------------------------------
IMPLICIT NONE
!-
CHARACTER(LEN=7) :: fun
INTEGER :: nb(2),nb_max,nb_index
INTEGER :: nindex(nb_index)
REAL :: work_in(nb(1),nb(2)),scal,miss_val
REAL :: work_out(nb_max)
!-
INTEGER :: ierr
!-
INTRINSIC SIN,COS,TAN,ASIN,ACOS,ATAN,EXP,LOG,SQRT,ABS
!---------------------------------------------------------------------
ierr = 0
!-
IF (scal >= miss_val-1.) THEN
IF (INDEX(indexfu,fun(1:LEN_TRIM(fun))) == 0) THEN
SELECT CASE (fun)
CASE('sin')
ierr = ma_sin_r21(nb,work_in,nb_max,work_out)
CASE('cos')
ierr = ma_cos_r21(nb,work_in,nb_max,work_out)
CASE('tan')
ierr = ma_tan_r21(nb,work_in,nb_max,work_out)
CASE('asin')
ierr = ma_asin_r21(nb,work_in,nb_max,work_out)
CASE('acos')
ierr = ma_acos_r21(nb,work_in,nb_max,work_out)
CASE('atan')
ierr = ma_atan_r21(nb,work_in,nb_max,work_out)
CASE('exp')
ierr = ma_exp_r21(nb,work_in,nb_max,work_out)
CASE('log')
ierr = ma_log_r21(nb,work_in,nb_max,work_out)
CASE('sqrt')
ierr = ma_sqrt_r21(nb,work_in,nb_max,work_out)
CASE('chs')
ierr = ma_chs_r21(nb,work_in,nb_max,work_out)
CASE('abs')
ierr = ma_abs_r21(nb,work_in,nb_max,work_out)
CASE('cels')
ierr = ma_cels_r21(nb,work_in,nb_max,work_out)
CASE('kelv')
ierr = ma_kelv_r21(nb,work_in,nb_max,work_out)
CASE('deg')
ierr = ma_deg_r21(nb,work_in,nb_max,work_out)
CASE('rad')
ierr = ma_rad_r21(nb,work_in,nb_max,work_out)
CASE('ident')
ierr = ma_ident_r21(nb,work_in,nb_max,work_out)
CASE DEFAULT
CALL ipslerr(3,"mathop", &
& 'scalar variable undefined and no indexing', &
& 'but still unknown function',fun)
END SELECT
IF (ierr > 0) THEN
CALL ipslerr(3,"mathop", &
& 'Error while executing a simple function',fun,' ')
ENDIF
ELSE
SELECT CASE (fun)
CASE('gather')
ierr = ma_fugath_r21(nb,work_in,nb_index,nindex, &
& miss_val,nb_max,work_out)
CASE('scatter')
IF (nb_index > (nb(1)*nb(2)) ) THEN
work_out(1:nb_max) = miss_val
ierr=1
ELSE
ierr = ma_fuscat_r21(nb,work_in,nb_index,nindex, &
& miss_val,nb_max,work_out)
ENDIF
CASE('coll')
ierr = ma_fucoll_r21(nb,work_in,nb_index,nindex, &
& miss_val,nb_max,work_out)
CASE('fill')
ierr = ma_fufill_r21(nb,work_in,nb_index,nindex, &
& miss_val,nb_max,work_out)
CASE('undef')
ierr = ma_fuundef_r21(nb,work_in,nb_index,nindex, &
& miss_val,nb_max,work_out)
CASE('only')
ierr = ma_fuonly_r21(nb,work_in,nb_index,nindex, &
& miss_val,nb_max,work_out)
CASE DEFAULT
CALL ipslerr(3,"mathop", &
& 'scalar variable undefined and indexing', &
& 'was requested but with unknown function',fun)
END SELECT
IF (ierr > 0) THEN
CALL ipslerr(3,"mathop_r21", &
& 'Error while executing an indexing function',fun,' ')
ENDIF
ENDIF
ELSE
SELECT CASE (fun)
CASE('fumin')
ierr = ma_fumin_r21(nb,work_in,scal,nb_max,work_out)
CASE('fumax')
ierr = ma_fumax_r21(nb,work_in,scal,nb_max,work_out)
CASE('add')
ierr = ma_add_r21(nb,work_in,scal,nb_max,work_out)
CASE('subi')
ierr = ma_subi_r21(nb,work_in,scal,nb_max,work_out)
CASE('sub')
ierr = ma_sub_r21(nb,work_in,scal,nb_max,work_out)
CASE('mult')
ierr = ma_mult_r21(nb,work_in,scal,nb_max,work_out)
CASE('div')
ierr = ma_div_r21(nb,work_in,scal,nb_max,work_out)
CASE('divi')
ierr = ma_divi_r21(nb,work_in,scal,nb_max,work_out)
CASE('power')
ierr = ma_power_r21(nb,work_in,scal,nb_max,work_out)
CASE DEFAULT
CALL ipslerr(3,"mathop", &
& 'Unknown operation with a scalar',fun,' ')
END SELECT
IF (ierr > 0) THEN
CALL ipslerr(3,"mathop", &
& 'Error while executing a scalar function',fun,' ')
ENDIF
ENDIF
!------------------------
END SUBROUTINE mathop_r21
!-
!=== FUNCTIONS (only one argument)
!-
INTEGER FUNCTION ma_sin_r21(nb,x,nbo,y)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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_log_r21(nb,x,nbo,y)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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) = LOG(x(i,j))
ENDDO
ENDDO
!-
nbo = nb(1)*nb(2)
ma_log_r21 = 0
!----------------------
END FUNCTION ma_log_r21
!===
INTEGER FUNCTION ma_sqrt_r21(nb,x,nbo,y)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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
!-----------------------
END FUNCTION ma_kelv_r21
!===
INTEGER FUNCTION ma_deg_r21(nb,x,nbo,y)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)*57.29577951
ENDDO
ENDDO
!-
nbo = nb(1)*nb(2)
ma_deg_r21 = 0
!----------------------
END FUNCTION ma_deg_r21
!===
INTEGER FUNCTION ma_rad_r21(nb,x,nbo,y)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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
!-
!- OUTPUT
!-
!- nb_max : Actual length of output variable
!- work_out : Output vector after the operation was applied
!---------------------------------------------------------------------
IMPLICIT NONE
!-
CHARACTER(LEN=7) :: fun
INTEGER :: nb(3),nb_max,nb_index
INTEGER :: nindex(nb_index)
REAL :: work_in(nb(1),nb(2),nb(3)),scal,miss_val
REAL :: work_out(nb_max)
!-
INTEGER :: ierr
!-
INTRINSIC SIN,COS,TAN,ASIN,ACOS,ATAN,EXP,LOG,SQRT,ABS
!---------------------------------------------------------------------
ierr = 0
!-
IF (scal >= miss_val-1.) THEN
IF (INDEX(indexfu,fun(1:LEN_TRIM(fun))) == 0) THEN
SELECT CASE (fun)
CASE('sin')
ierr = ma_sin_r31(nb,work_in,nb_max,work_out)
CASE('cos')
ierr = ma_cos_r31(nb,work_in,nb_max,work_out)
CASE('tan')
ierr = ma_tan_r31(nb,work_in,nb_max,work_out)
CASE('asin')
ierr = ma_asin_r31(nb,work_in,nb_max,work_out)
CASE('acos')
ierr = ma_acos_r31(nb,work_in,nb_max,work_out)
CASE('atan')
ierr = ma_atan_r31(nb,work_in,nb_max,work_out)
CASE('exp')
ierr = ma_exp_r31(nb,work_in,nb_max,work_out)
CASE('log')
ierr = ma_log_r31(nb,work_in,nb_max,work_out)
CASE('sqrt')
ierr = ma_sqrt_r31(nb,work_in,nb_max,work_out)
CASE('chs')
ierr = ma_chs_r31(nb,work_in,nb_max,work_out)
CASE('abs')
ierr = ma_abs_r31(nb,work_in,nb_max,work_out)
CASE('cels')
ierr = ma_cels_r31(nb,work_in,nb_max,work_out)
CASE('kelv')
ierr = ma_kelv_r31(nb,work_in,nb_max,work_out)
CASE('deg')
ierr = ma_deg_r31(nb,work_in,nb_max,work_out)
CASE('rad')
ierr = ma_rad_r31(nb,work_in,nb_max,work_out)
CASE('ident')
ierr = ma_ident_r31(nb,work_in,nb_max,work_out)
CASE DEFAULT
CALL ipslerr(3,"mathop", &
& 'scalar variable undefined and no indexing', &
& 'but still unknown function',fun)
END SELECT
IF (ierr > 0) THEN
CALL ipslerr(3,"mathop", &
& 'Error while executing a simple function',fun,' ')
ENDIF
ELSE
SELECT CASE (fun)
CASE('gather')
ierr = ma_fugath_r31(nb,work_in,nb_index,nindex, &
& miss_val,nb_max,work_out)
CASE('scatter')
IF (nb_index > (nb(1)*nb(2)*nb(3))) THEN
work_out(1:nb_max) = miss_val
ierr=1
ELSE
ierr = ma_fuscat_r31(nb,work_in,nb_index,nindex, &
& miss_val,nb_max,work_out)
ENDIF
CASE('coll')
ierr = ma_fucoll_r31(nb,work_in,nb_index,nindex, &
& miss_val,nb_max,work_out)
CASE('fill')
ierr = ma_fufill_r31(nb,work_in,nb_index,nindex, &
& miss_val,nb_max,work_out)
CASE('undef')
ierr = ma_fuundef_r31(nb,work_in,nb_index,nindex, &
& miss_val,nb_max,work_out)
CASE('only')
ierr = ma_fuonly_r31(nb,work_in,nb_index,nindex, &
& miss_val,nb_max,work_out)
CASE DEFAULT
CALL ipslerr(3,"mathop", &
& 'scalar variable undefined and indexing', &
& 'was requested but with unknown function',fun)
END SELECT
IF (ierr > 0) THEN
CALL ipslerr(3,"mathop_r31", &
& 'Error while executing an indexing function',fun,' ')
ENDIF
ENDIF
ELSE
SELECT CASE (fun)
CASE('fumin')
ierr = ma_fumin_r31(nb,work_in,scal,nb_max,work_out)
CASE('fumax')
ierr = ma_fumax_r31(nb,work_in,scal,nb_max,work_out)
CASE('add')
ierr = ma_add_r31(nb,work_in,scal,nb_max,work_out)
CASE('subi')
ierr = ma_subi_r31(nb,work_in,scal,nb_max,work_out)
CASE('sub')
ierr = ma_sub_r31(nb,work_in,scal,nb_max,work_out)
CASE('mult')
ierr = ma_mult_r31(nb,work_in,scal,nb_max,work_out)
CASE('div')
ierr = ma_div_r31(nb,work_in,scal,nb_max,work_out)
CASE('divi')
ierr = ma_divi_r31(nb,work_in,scal,nb_max,work_out)
CASE('power')
ierr = ma_power_r31(nb,work_in,scal,nb_max,work_out)
CASE DEFAULT
CALL ipslerr(3,"mathop", &
& 'Unknown operation with a scalar',fun,' ')
END SELECT
IF (ierr > 0) THEN
CALL ipslerr(3,"mathop", &
& 'Error while executing a scalar function',fun,' ')
ENDIF
ENDIF
!------------------------
END SUBROUTINE mathop_r31
!-
!=== FUNCTIONS (only one argument)
!-
INTEGER FUNCTION ma_sin_r31(nb,x,nbo,y)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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
!===
INTEGER FUNCTION ma_exp_r31(nb,x,nbo,y)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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_log_r31(nb,x,nbo,y)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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) = LOG(x(i,j,k))
ENDDO
ENDDO
ENDDO
!-
nbo = nb(1)*nb(2)*nb(3)
ma_log_r31 = 0
!----------------------
END FUNCTION ma_log_r31
!===
INTEGER FUNCTION ma_sqrt_r31(nb,x,nbo,y)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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)
!---------------------------------------------------------------------
IMPLICIT NONE
!-
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
!===
SUBROUTINE moycum (opp,np,px,py,pwx)
!---------------------------------------------------------------------
!- Does time operations
!---------------------------------------------------------------------
IMPLICIT NONE
!-
CHARACTER(LEN=7) :: opp
INTEGER :: np
REAL,DIMENSION(:) :: px,py
INTEGER :: pwx
!---------------------------------------------------------------------
IF (pwx /= 0) THEN
IF (opp == 'ave') THEN
px(1:np)=(px(1:np)*pwx+py(1:np))/REAL(pwx+1)
ELSE IF (opp == 't_sum') THEN
px(1:np)=px(1:np)+py(1:np)
ELSE IF ( (opp == 'l_min').OR.(opp == 't_min') ) THEN
px(1:np)=MIN(px(1:np),py(1:np))
ELSE IF ( (opp == 'l_max').OR.(opp == 't_max') ) THEN
px(1:np)=MAX(px(1:np),py(1:np))
ELSE
CALL ipslerr(3,"moycum",'Unknown time operation',opp,' ')
ENDIF
ELSE
IF (opp == 'l_min') THEN
px(1:np)=MIN(px(1:np),py(1:np))
ELSE IF (opp == 'l_max') THEN
px(1:np)=MAX(px(1:np),py(1:np))
ELSE
px(1:np)=py(1:np)
ENDIF
ENDIF
!--------------------
END SUBROUTINE moycum
!===
!-----------------
END MODULE mathelp
Reference in New Issue View Git Blame Copy Permalink