/[lmdze]/trunk/IOIPSL/Mathelp/mathelp.f
ViewVC logotype

Annotation of /trunk/IOIPSL/Mathelp/mathelp.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 32 - (hide annotations)
Tue Apr 6 17:52:58 2010 UTC (14 years, 2 months ago) by guez
Original Path: trunk/libf/IOIPSL/mathelp.f90
File size: 18105 byte(s)
Split "stringop.f90" into single-procedure files. Gathered files in directory
"IOIPSL/Stringop".

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

Removed unused files.

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

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

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

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

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

Removed unused arguments of procedure "integrd".

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

Merged common block "matrfil" into module "parafilt".

1 guez 30 MODULE mathelp
2    
3     ! From mathelp.f90, version 2.0 2004/04/05 14:47:50
4    
5     implicit none
6    
7     PRIVATE
8 guez 32 PUBLIC :: moycum, trans_buff, buildop
9 guez 30
10     !- Variables used to detect and identify the operations
11 guez 32 CHARACTER(LEN=80), SAVE :: seps='( ) , + - / * ^'
12 guez 30
13     CONTAINS
14    
15     SUBROUTINE buildop (str, ex_topps, topp, nbops_max, &
16     & missing_val, opps, scal, nbops)
17     !- This subroutine decomposes the input string in the elementary
18     !- functions which need to be applied to the vector of data.
19     !- This vector is represented by X in the string.
20     !- This subroutine is the driver of the decomposition and gets
21     !- the time operation but then call decoop for the other operations
22     !- INPUT
23    
24     !- str : String containing the operations
25     !- ex_toops : The time operations that can be expected
26     !- within the string
27    
28     !- OUTPUT
29    
30 guez 32 USE errioipsl, ONLY : histerr
31    
32 guez 30 CHARACTER(LEN=80) :: str
33     CHARACTER(LEN=*) :: ex_topps
34     CHARACTER(LEN=7) :: topp
35     INTEGER :: nbops_max, nbops
36     CHARACTER(LEN=7) :: opps(nbops_max)
37     REAL :: scal(nbops_max), missing_val
38    
39     CHARACTER(LEN=80) :: new_str
40     INTEGER :: leng, ind_opb, ind_clb
41    
42     LOGICAL :: check = .FALSE.
43     !---------------------------------------------------------------------
44     IF (check) WRITE(*, *) 'buildop : Some preliminary cleaning'
45    
46     leng = LEN_TRIM(str)
47     IF ( str(1:1) == '(' .AND. str(leng:leng) == ')' ) THEN
48     str = str(2:leng-1)
49     leng = leng-2
50     ENDIF
51    
52     IF (check) &
53     & WRITE(*, *) 'buildop : Starting to test the various options'
54    
55     IF (leng <= 5 .AND. INDEX(ex_topps, str(1:leng)) > 0) THEN
56     IF (check) WRITE(*, *) 'buildop : Time operation only'
57     nbops = 0
58     topp = str(1:leng)
59     ELSE
60     IF (check) THEN
61     WRITE(*, *) 'buildop : Time operation and something else'
62     ENDIF
63     !--
64     ind_opb = INDEX(str(1:leng), '(')
65     IF (ind_opb > 0) THEN
66     IF (INDEX(ex_topps, str(1:ind_opb-1)) > 0) THEN
67     IF (check) THEN
68     WRITE(*, '(2a)') &
69     & ' buildop : Extract time operation from : ', str
70     ENDIF
71     topp = str(1:ind_opb-1)
72     ind_clb = INDEX(str(1:leng), ')', BACK=.TRUE.)
73     new_str = str(ind_opb+1:ind_clb-1)
74     IF (check) THEN
75     WRITE(*, '(2a, 2I3)') &
76     & ' buildop : Call decoop ', new_str, ind_opb, ind_clb
77     ENDIF
78     CALL decoop (new_str, nbops_max, missing_val, opps, scal, nbops)
79     ELSE
80     CALL histerr(3, 'buildop', &
81     & 'time opperation does not exist', str(1:ind_opb-1), ' ')
82     ENDIF
83     ELSE
84     CALL histerr(3, 'buildop', &
85     & 'some long opperation exists but wihout parenthesis', &
86     & str(1:leng), ' ')
87     ENDIF
88     ENDIF
89    
90     IF (check) THEN
91     DO leng=1, nbops
92     WRITE(*, *) &
93     & 'buildop : i -- opps, scal : ', leng, opps(leng), scal(leng)
94     ENDDO
95     ENDIF
96     !---------------------
97     END SUBROUTINE buildop
98    
99     !************************************************
100    
101     SUBROUTINE decoop (pstr, nbops_max, missing_val, opps, scal, nbops)
102 guez 32 USE errioipsl, ONLY : histerr
103    
104 guez 30 CHARACTER(LEN=80) :: pstr
105     INTEGER :: nbops_max, nbops
106     CHARACTER(LEN=7) :: opps(nbops_max)
107     REAL :: scal(nbops_max), missing_val
108    
109     CHARACTER(LEN=1) :: f_char(2), s_char(2)
110     INTEGER :: nbsep, f_pos(2), s_pos(2)
111     CHARACTER(LEN=20) :: opp_str, scal_str
112     CHARACTER(LEN=80) :: str
113     INTEGER :: xpos, leng, ppos, epos, int_tmp
114     CHARACTER(LEN=3) :: tl, dl
115     CHARACTER(LEN=10) :: fmt
116    
117     LOGICAL :: check = .FALSE., prio
118 guez 32 CHARACTER(LEN=80), SAVE :: ops = '+ - * / ^'
119     CHARACTER(LEN=80), SAVE :: mima = 'min max'
120     CHARACTER(LEN=250), SAVE :: funcs = &
121     'sin cos tan asin acos atan exp log sqrt chs abs ' &
122     //'cels kelv deg rad gather scatter fill coll undef only ident'
123 guez 30 !---------------------------------------------------------------------
124     IF (check) WRITE(*, '(2a)') ' decoop : Incoming string : ', pstr
125    
126     nbops = 0
127     str = pstr
128    
129     CALL findsep (str, nbsep, f_char, f_pos, s_char, s_pos)
130     IF (check) WRITE(*, *) 'decoop : Out of findsep', nbsep
131     DO WHILE (nbsep > 0)
132     xpos = INDEX(str, 'X')
133     leng = LEN_TRIM(str)
134     nbops = nbops+1
135     !--
136     IF (check) THEN
137     WRITE(*, *) 'decoop : str -->', str(1:leng)
138     WRITE(*, *) s_char(1), '-', f_char(1), '|', f_char(2), '-', s_char(2)
139     WRITE(*, *) s_pos(1), '-', f_pos(1), '|', f_pos(2), '-', s_pos(2)
140     ENDIF
141     !--
142     IF (nbops > nbops_max-1) THEN
143     CALL histerr(3, 'decoop', 'Expression too complex', str, ' ')
144     ENDIF
145     !--
146     IF (check) WRITE(*, *) 'decoop : --', nbops, ' ', str(1:leng)
147     !---
148     !-- Start the analysis of the syntax. 3 types of constructs
149     !-- are recognized. They are scanned sequentialy
150     !---
151     IF (nbsep == 1) THEN
152     IF (check) WRITE(*, *) 'decoop : Only one operation'
153     IF (INDEX(ops, f_char(1)) > 0) THEN
154     !------ Type : scal+X
155     IF (f_char(1) == '-' .OR. f_char(1) == '/') THEN
156     opp_str = f_char(1)//'I'
157     ELSE
158     opp_str = f_char(1)
159     ENDIF
160     scal_str = str(s_pos(1)+1:f_pos(1)-1)
161     str = 'X'
162     ELSE IF (INDEX(ops, f_char(2)) > 0) THEN
163     !------ Type : X+scal
164     opp_str = f_char(2)
165     scal_str = str(f_pos(2)+1:s_pos(2)-1)
166     str = 'X'
167     ELSE
168     CALL histerr(3, 'decoop', &
169     & 'Unknown operations of type X+scal', f_char(1), pstr)
170     ENDIF
171     ELSE
172     IF (check) WRITE(*, *) 'decoop : More complex operation'
173     IF ( f_char(1) == '(' .AND. f_char(2) == ')' ) THEN
174     !------ Type : sin(X)
175     opp_str = str(s_pos(1)+1:f_pos(1)-1)
176     scal_str = '?'
177     str = str(1:s_pos(1))//'X'//str(f_pos(2)+1:leng)
178     ELSE IF ( (f_char(1) == '(' .AND. f_char(2) == ', ')&
179     & .OR.(f_char(1) == ', ' .AND. f_char(2) == ')')) THEN
180     !------ Type : max(X, scal) or max(scal, X)
181     IF (f_char(1) == '(' .AND. s_char(2) == ')') THEN
182     !-------- Type : max(X, scal)
183     opp_str = str(f_pos(1)-3:f_pos(1)-1)
184     scal_str = str(f_pos(2)+1:s_pos(2)-1)
185     str = str(1:f_pos(1)-4)//'X'//str(s_pos(2)+1:leng)
186     ELSE IF (f_char(1) == ', ' .AND. s_char(1) == '(') THEN
187     !-------- Type : max(scal, X)
188     opp_str = str(s_pos(1)-3:s_pos(1)-1)
189     scal_str = str(s_pos(1)+1:f_pos(1)-1)
190     str = str(1:s_pos(1)-4)//'X'//str(f_pos(2)+1:leng)
191     ELSE
192     CALL histerr(3, 'decoop', 'Syntax error 1', str, ' ')
193     ENDIF
194     ELSE
195     prio = (f_char(2) == '*').OR.(f_char(2) == '^')
196     IF ( (INDEX(ops, f_char(1)) > 0) &
197     & .AND.(xpos-f_pos(1) == 1).AND.(.NOT.prio) ) THEN
198     !-------- Type : ... scal+X ...
199     IF (f_char(1) == '-' .OR. f_char(1) == '/') THEN
200     opp_str = f_char(1)//'I'
201     ELSE
202     opp_str = f_char(1)
203     ENDIF
204     scal_str = str(s_pos(1)+1:f_pos(1)-1)
205     str = str(1:s_pos(1))//'X'//str(f_pos(1)+2:leng)
206     ELSE IF ( (INDEX(ops, f_char(2)) > 0) &
207     & .AND.(f_pos(2)-xpos == 1) ) THEN
208     !-------- Type : ... X+scal ...
209     opp_str = f_char(2)
210     scal_str = str(f_pos(2)+1:s_pos(2)-1)
211     str = str(1:f_pos(2)-2)//'X'//str(s_pos(2):leng)
212     ELSE
213     CALL histerr(3, 'decoop', 'Syntax error 2', str, ' ')
214     ENDIF
215     ENDIF
216     ENDIF
217     !---
218     IF (check) WRITE(*, *) 'decoop : Finished syntax, str = ', TRIM(str)
219     !---
220     !-- Now that the different components of the operation are identified
221     !-- we transform them into what is going to be used in the program
222     !---
223     IF (INDEX(scal_str, '?') > 0) THEN
224     IF (INDEX(funcs, opp_str(1:LEN_TRIM(opp_str))) > 0) THEN
225     opps(nbops) = opp_str(1:LEN_TRIM(opp_str))
226     scal(nbops) = missing_val
227     ELSE
228     CALL histerr(3, 'decoop', &
229     & 'Unknown function', opp_str(1:LEN_TRIM(opp_str)), ' ')
230     ENDIF
231     ELSE
232     leng = LEN_TRIM(opp_str)
233     IF (INDEX(mima, opp_str(1:leng)) > 0) THEN
234     opps(nbops) = 'fu'//opp_str(1:leng)
235     ELSE
236     IF (INDEX(opp_str(1:leng), '+') > 0) THEN
237     opps(nbops) = 'add'
238     ELSE IF (INDEX(opp_str(1:leng), '-I') > 0) THEN
239     opps(nbops) = 'subi'
240     ELSE IF (INDEX(opp_str(1:leng), '-') > 0) THEN
241     opps(nbops) = 'sub'
242     ELSE IF (INDEX(opp_str(1:leng), '*') > 0) THEN
243     opps(nbops) = 'mult'
244     ELSE IF (INDEX(opp_str(1:leng), '/') > 0) THEN
245     opps(nbops) = 'div'
246     ELSE IF (INDEX(opp_str(1:leng), '/I') > 0) THEN
247     opps(nbops) = 'divi'
248     ELSE IF (INDEX(opp_str(1:leng), '^') > 0) THEN
249     opps(nbops) = 'power'
250     ELSE
251     CALL histerr(3, 'decoop', &
252     & 'Unknown operation', opp_str(1:leng), ' ')
253     ENDIF
254     ENDIF
255     !-----
256     leng = LEN_TRIM(scal_str)
257     ppos = INDEX(scal_str, '.')
258     epos = INDEX(scal_str, 'e')
259     IF (epos == 0) epos = INDEX(scal_str, 'E')
260     !-----
261     !---- Try to catch a few errors
262     !-----
263     IF (INDEX(ops, scal_str) > 0) THEN
264     CALL histerr(3, 'decoop', &
265     & 'Strange scalar you have here ', scal_str, pstr)
266     ENDIF
267     IF (epos > 0) THEN
268     WRITE(tl, '(I3.3)') leng
269     WRITE(dl, '(I3.3)') epos-ppos-1
270     fmt='(e'//tl//'.'//dl//')'
271     READ(scal_str, fmt) scal(nbops)
272     ELSE IF (ppos > 0) THEN
273     WRITE(tl, '(I3.3)') leng
274     WRITE(dl, '(I3.3)') leng-ppos
275     fmt='(f'//tl//'.'//dl//')'
276     READ(scal_str, fmt) scal(nbops)
277     ELSE
278     WRITE(tl, '(I3.3)') leng
279     fmt = '(I'//tl//')'
280     READ(scal_str, fmt) int_tmp
281     scal(nbops) = REAL(int_tmp)
282     ENDIF
283     ENDIF
284     IF (check) WRITE(*, *) 'decoop : Finished interpretation'
285     CALL findsep(str, nbsep, f_char, f_pos, s_char, s_pos)
286     ENDDO
287     !--------------------
288     END SUBROUTINE decoop
289    
290     !************************************************
291    
292     SUBROUTINE findsep (str, nbsep, f_char, f_pos, s_char, s_pos)
293     !- Subroutine finds all separators in a given string
294     !- It returns the following information about str :
295     !- f_char : The first separation character
296     !- (1 for before and 2 for after)
297     !- f_pos : The position of the first separator
298     !- s_char : The second separation character
299     !- (1 for before and 2 for after)
300     !- s_pos : The position of the second separator
301 guez 32 USE errioipsl, ONLY : histerr
302    
303 guez 30 CHARACTER(LEN=80) :: str
304     INTEGER :: nbsep
305     CHARACTER(LEN=1), DIMENSION(2) :: f_char, s_char
306     INTEGER, DIMENSION(2) :: f_pos, s_pos
307    
308     CHARACTER(LEN=70) :: str_tmp
309     LOGICAL :: f_found, s_found
310     INTEGER :: ind, xpos, leng, i
311    
312     LOGICAL :: check = .FALSE.
313     !---------------------------------------------------------------------
314     IF (check) WRITE(*, *) 'findsep : call cleanstr: ', TRIM(str)
315    
316     CALL cleanstr(str)
317    
318     IF (check) WRITE(*, *) 'findsep : out of cleanstr: ', TRIM(str)
319    
320     xpos = INDEX(str, 'X')
321     leng = LEN_TRIM(str)
322    
323     f_pos(1:2) = (/ 0, leng+1 /)
324     f_char(1:2) = (/ '?', '?' /)
325     s_pos(1:2) = (/ 0, leng+1 /)
326     s_char(1:2) = (/ '?', '?' /)
327    
328     nbsep = 0
329    
330     f_found = .FALSE.
331     s_found = .FALSE.
332     IF (xpos > 1) THEN
333     DO i=xpos-1, 1, -1
334     ind = INDEX(seps, str(i:i))
335     IF (ind > 0) THEN
336     IF (.NOT.f_found) THEN
337     f_char(1) = str(i:i)
338     f_pos(1) = i
339     nbsep = nbsep+1
340     f_found = .TRUE.
341     ELSE IF (.NOT.s_found) THEN
342     s_char(1) = str(i:i)
343     s_pos(1) = i
344     nbsep = nbsep+1
345     s_found = .TRUE.
346     ENDIF
347     ENDIF
348     ENDDO
349     ENDIF
350    
351     f_found = .FALSE.
352     s_found = .FALSE.
353     IF (xpos < leng) THEN
354     DO i=xpos+1, leng
355     ind = INDEX(seps, str(i:i))
356     IF (ind > 0) THEN
357     IF (.NOT.f_found) THEN
358     f_char(2) = str(i:i)
359     f_pos(2) = i
360     nbsep = nbsep+1
361     f_found = .TRUE.
362     ELSE IF (.NOT.s_found) THEN
363     s_char(2) = str(i:i)
364     s_pos(2) = i
365     nbsep = nbsep+1
366     s_found = .TRUE.
367     ENDIF
368     ENDIF
369     ENDDO
370     ENDIF
371    
372     IF (nbsep > 4) THEN
373     WRITE(str_tmp, '("number :", I3)') nbsep
374     CALL histerr(3, 'findsep', &
375     & 'How can I find that many separators', str_tmp, str)
376     ENDIF
377    
378     IF (check) WRITE(*, *) 'Finished findsep : ', nbsep, leng
379     !---------------------
380     END SUBROUTINE findsep
381    
382     !************************************************
383    
384     SUBROUTINE cleanstr(str)
385     !- We clean up the string by taking out the extra () and puting
386     !- everything in lower case except for the X describing the variable
387 guez 32 use strlowercase_m, only: strlowercase
388    
389 guez 30 CHARACTER(LEN=80) :: str
390    
391     INTEGER :: ind, leng, ic, it
392     LOGICAL :: check = .FALSE.
393     !---------------------------------------------------------------------
394     leng = LEN_TRIM(str)
395     CALL strlowercase(str)
396    
397     ind = INDEX(str, 'x')
398     IF (check) THEN
399     WRITE (*, *) 'cleanstr 1.0 : ind = ', ind, &
400     & ' str = ', str(1:leng), '---'
401     ENDIF
402    
403     ! If the character before the x is not a letter then we can assume
404     ! that it is the variable and promote it to a capital letter
405    
406     DO WHILE (ind > 0)
407     ic = 0
408     IF (ind > 1) ic = IACHAR(str(ind-1:ind-1))
409     IF (ic < 97 .OR. ic > 122) THEN
410     str(ind:ind) = 'X'
411     ENDIF
412     it = INDEX(str(ind+1:leng), 'x')
413     IF (it > 0) THEN
414     ind = ind+it
415     ELSE
416     ind = it
417     ENDIF
418     ENDDO
419    
420     IF (check) WRITE (*, *) 'cleanstr 2.0 : str = ', str(1:leng), '---'
421    
422     IF ( str(1:1) == '(' .AND. str(leng:leng) == ')' ) THEN
423     str = str(2:leng-1)
424     ENDIF
425    
426     IF (check) WRITE (*, *) 'cleanstr 3.0 : str = ', str(1:leng), '---'
427    
428     leng = LEN_TRIM(str)
429     ind = INDEX(str, '((X))')
430     IF (ind > 0) THEN
431     str=str(1:ind-1)//'(X)'//str(ind+5:leng)//' '
432     ENDIF
433    
434     IF (check) WRITE (*, *) 'cleanstr 4.0 : str = ', str(1:leng), '---'
435    
436     leng = LEN_TRIM(str)
437     ind = INDEX(str, '(X)')
438     IF (ind > 0 .AND. ind+3 < leng) THEN
439     IF ( (INDEX(seps, str(ind-1:ind-1)) > 0) &
440     & .AND. (INDEX(seps, str(ind+3:ind+3)) > 0) ) THEN
441     str=str(1:ind-1)//'X'//str(ind+3:leng)//' '
442     ENDIF
443     ENDIF
444    
445     IF (check) WRITE (*, *) 'cleanstr 5.0 : str = ', str(1:leng), '---'
446    
447     leng = LEN_TRIM(str)
448     ind = INDEX(str(1:leng), ' ')
449     DO WHILE (ind > 0)
450     str=str(1:ind-1)//str(ind+1:leng)//' '
451     leng = LEN_TRIM(str)
452     ind = INDEX(str(1:leng), ' ')
453     ENDDO
454    
455     IF (check) WRITE (*, *) 'cleanstr 6.0 : str = ', str(1:leng), '---'
456     !----------------------
457     END SUBROUTINE cleanstr
458    
459     !************************************************
460    
461     SUBROUTINE moycum (opp, np, px, py, pwx)
462     !- Does time operations
463 guez 32 USE errioipsl, ONLY : histerr
464    
465 guez 30 CHARACTER(LEN=7) :: opp
466     INTEGER :: np
467     REAL, DIMENSION(:) :: px, py
468     INTEGER :: pwx
469     !---------------------------------------------------------------------
470     IF (pwx /= 0) THEN
471     IF (opp == 'ave') THEN
472     px(1:np)=(px(1:np)*pwx+py(1:np))/REAL(pwx+1)
473     ELSE IF (opp == 't_sum') THEN
474     px(1:np)=px(1:np)+py(1:np)
475     ELSE IF ( (opp == 'l_min').OR.(opp == 't_min') ) THEN
476     px(1:np)=MIN(px(1:np), py(1:np))
477     ELSE IF ( (opp == 'l_max').OR.(opp == 't_max') ) THEN
478     px(1:np)=MAX(px(1:np), py(1:np))
479     ELSE
480     CALL histerr(3, "moycum", 'Unknown time operation', opp, ' ')
481     ENDIF
482     ELSE
483     IF (opp == 'l_min') THEN
484     px(1:np)=MIN(px(1:np), py(1:np))
485     ELSE IF (opp == 'l_max') THEN
486     px(1:np)=MAX(px(1:np), py(1:np))
487     ELSE
488     px(1:np)=py(1:np)
489     ENDIF
490     ENDIF
491     !--------------------
492     END SUBROUTINE moycum
493    
494     !************************************************
495    
496     SUBROUTINE trans_buff (ox, sx, oy, sy, oz, sz, xsz, ysz, zsz, v3d, sl, v1d)
497     !- This subroutine extracts from the full 3D variable the slab of
498     !- data that will be used later. Perhaps there are hardware routines
499     !- for this task on some computers. This routine will be obsolete in
500     !- a F90 environnement
501    
502     !- INPUT
503     !- ox : Origin of slab of data in X
504     !- sx : Size of slab in X
505     !- oy : Origin of slab of data in Y
506     !- sy : Size of slab in Y
507     !- oz : Origin of slab of data in Z
508     !- sz : Size of slab in Z
509     !- xsz, ysz, zsz : 3 sizes of full variable v3d
510     !- v3d : The full 3D variable
511     !- sl : size of variable v1d
512     !- v1d : The 1D variable containing the slab
513    
514     INTEGER :: ox, sx, oy, sy, oz, sz
515     INTEGER :: xsz, ysz, zsz
516     INTEGER :: sl
517     REAL :: v3d(xsz, ysz, zsz)
518     REAL :: v1d(sl)
519    
520     !---------------------------------------------------------------------
521    
522     V1d(:sx*sy*sz) = reshape(V3d(ox:ox-1+sx, oy:oy-1+sy, oz:oz-1+sz), &
523     SHAPE=(/sx*sy*sz/))
524    
525     END SUBROUTINE trans_buff
526    
527     END MODULE mathelp

  ViewVC Help
Powered by ViewVC 1.1.21