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

Contents of /trunk/libf/IOIPSL/mathelp.f90

Parent Directory Parent Directory | Revision Log Revision Log


Revision 32 - (show annotations)
Tue Apr 6 17:52:58 2010 UTC (14 years, 2 months ago) by guez
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 MODULE mathelp
2
3 ! From mathelp.f90, version 2.0 2004/04/05 14:47:50
4
5 implicit none
6
7 PRIVATE
8 PUBLIC :: moycum, trans_buff, buildop
9
10 !- Variables used to detect and identify the operations
11 CHARACTER(LEN=80), SAVE :: seps='( ) , + - / * ^'
12
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 USE errioipsl, ONLY : histerr
31
32 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 USE errioipsl, ONLY : histerr
103
104 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 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 !---------------------------------------------------------------------
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 USE errioipsl, ONLY : histerr
302
303 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 use strlowercase_m, only: strlowercase
388
389 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 USE errioipsl, ONLY : histerr
464
465 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