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

Contents of /trunk/IOIPSL/mathelp.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 82 - (show annotations)
Wed Mar 5 14:57:53 2014 UTC (10 years, 2 months ago) by guez
File size: 18105 byte(s)
Changed all ".f90" suffixes to ".f".
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