My Project
 All Classes Files Functions Variables Macros
grid_atob.F
Go to the documentation of this file.
1 !
2 ! $Id: grid_atob.F 1403 2010-07-01 09:02:53Z fairhead $
3 !
4  SUBROUTINE grille_m(imdep, jmdep, xdata, ydata, entree,
5  . imar, jmar, x, y, sortie)
6 c=======================================================================
7 c z.x.li (le 1 avril 1994) (voir aussi A. Harzallah et L. Fairhead)
8 c
9 c Methode naive pour transformer un champ d'une grille fine a une
10 c grille grossiere. Je considere que les nouveaux points occupent
11 c une zone adjacente qui comprend un ou plusieurs anciens points
12 c
13 c Aucune ponderation est consideree (voir grille_p)
14 c
15 c (c)
16 c ----d-----
17 c | . . . .|
18 c | |
19 c (b)a . * . .b(a)
20 c | |
21 c | . . . .|
22 c ----c-----
23 c (d)
24 C=======================================================================
25 c INPUT:
26 c imdep, jmdep: dimensions X et Y pour depart
27 c xdata, ydata: coordonnees X et Y pour depart
28 c entree: champ d'entree a transformer
29 c OUTPUT:
30 c imar, jmar: dimensions X et Y d'arrivee
31 c x, y: coordonnees X et Y d'arrivee
32 c sortie: champ de sortie deja transforme
33 C=======================================================================
34  IMPLICIT none
35 
36  INTEGER imdep, jmdep
37  REAL xdata(imdep),ydata(jmdep)
38  REAL entree(imdep,jmdep)
39 c
40  INTEGER imar, jmar
41  REAL x(imar),y(jmar)
42  REAL sortie(imar,jmar)
43 c
44  INTEGER i, j, ii, jj
45  REAL a(2200),b(2200),c(1100),d(1100)
46  REAL number(2200,1100)
47  REAL distans(2200*1100)
48  INTEGER i_proche, j_proche, ij_proche
49 #ifdef CRAY
50  INTEGER ismin
51 #else
52  REAL zzmin
53 #endif
54 c
55  IF (imar.GT.2200 .OR. jmar.GT.1100) THEN
56  print*, 'imar ou jmar trop grand', imar, jmar
57  CALL abort
58  ENDIF
59 c
60 c Calculer les limites des zones des nouveaux points
61 c
62 
63  a(1) = x(1) - (x(2)-x(1))/2.0
64  b(1) = (x(1)+x(2))/2.0
65  DO i = 2, imar-1
66  a(i) = b(i-1)
67  b(i) = (x(i)+x(i+1))/2.0
68  ENDDO
69  a(imar) = b(imar-1)
70  b(imar) = x(imar) + (x(imar)-x(imar-1))/2.0
71 
72  c(1) = y(1) - (y(2)-y(1))/2.0
73  d(1) = (y(1)+y(2))/2.0
74  DO j = 2, jmar-1
75  c(j) = d(j-1)
76  d(j) = (y(j)+y(j+1))/2.0
77  ENDDO
78  c(jmar) = d(jmar-1)
79  d(jmar) = y(jmar) + (y(jmar)-y(jmar-1))/2.0
80 
81  DO i = 1, imar
82  DO j = 1, jmar
83  number(i,j) = 0.0
84  sortie(i,j) = 0.0
85  ENDDO
86  ENDDO
87 c
88 c Determiner la zone sur laquelle chaque ancien point se trouve
89 c
90 c
91 c ..... Modif P. Le Van ( 23/08/95 ) ....
92 
93  DO ii = 1, imar
94  DO jj = 1, jmar
95  DO i = 1, imdep
96  IF( ( xdata(i)-a(ii).GE.1.e-5.AND.xdata(i)-b(ii).LE.1.e-5 ).OR.
97  . ( xdata(i)-a(ii).LE.1.e-5.AND.xdata(i)-b(ii).GE.1.e-5 ) )
98  . THEN
99  DO j = 1, jmdep
100  IF( (ydata(j)-c(jj).GE.1.e-5.AND.ydata(j)-d(jj).LE.1.e-5 ).OR.
101  . ( ydata(j)-c(jj).LE.1.e-5.AND.ydata(j)-d(jj).GE.1.e-5 ) )
102  . THEN
103  number(ii,jj) = number(ii,jj) + 1.0
104  sortie(ii,jj) = sortie(ii,jj) + entree(i,j)
105  ENDIF
106  ENDDO
107  ENDIF
108  ENDDO
109  ENDDO
110  ENDDO
111 c
112 c Si aucun ancien point tombe sur une zone, c'est un probleme
113 c
114 
115  DO i = 1, imar
116  DO j = 1, jmar
117  IF (number(i,j) .GT. 0.001) THEN
118  sortie(i,j) = sortie(i,j) / number(i,j)
119  ELSE
120  print*, 'probleme,i,j=', i,j
121 ccc CALL ABORT
122  CALL dist_sphe(x(i),y(j),xdata,ydata,imdep,jmdep,distans)
123 #ifdef CRAY
124  ij_proche = ismin(imdep*jmdep,distans,1)
125 #else
126  ij_proche = 1
127  zzmin = distans(ij_proche)
128  DO ii = 2, imdep*jmdep
129  IF (distans(ii).LT.zzmin) THEN
130  zzmin = distans(ii)
131  ij_proche = ii
132  ENDIF
133  ENDDO
134 #endif
135  j_proche = (ij_proche-1)/imdep + 1
136  i_proche = ij_proche - (j_proche-1)*imdep
137  print*, "solution:", ij_proche, i_proche, j_proche
138  sortie(i,j) = entree(i_proche,j_proche)
139  ENDIF
140  ENDDO
141  ENDDO
142 
143  RETURN
144  END
145 
146 
147  SUBROUTINE grille_p(imdep, jmdep, xdata, ydata, entree,
148  . imar, jmar, x, y, sortie)
149 c=======================================================================
150 c z.x.li (le 1 avril 1994) (voir aussi A. Harzallah et L. Fairhead)
151 c
152 c Methode naive pour transformer un champ d'une grille fine a une
153 c grille grossiere. Je considere que les nouveaux points occupent
154 c une zone adjacente qui comprend un ou plusieurs anciens points
155 c
156 c Consideration de la distance des points (voir grille_m)
157 c
158 c (c)
159 c ----d-----
160 c | . . . .|
161 c | |
162 c (b)a . * . .b(a)
163 c | |
164 c | . . . .|
165 c ----c-----
166 c (d)
167 C=======================================================================
168 c INPUT:
169 c imdep, jmdep: dimensions X et Y pour depart
170 c xdata, ydata: coordonnees X et Y pour depart
171 c entree: champ d'entree a transformer
172 c OUTPUT:
173 c imar, jmar: dimensions X et Y d'arrivee
174 c x, y: coordonnees X et Y d'arrivee
175 c sortie: champ de sortie deja transforme
176 C=======================================================================
177  IMPLICIT none
178 
179  INTEGER imdep, jmdep
180  REAL xdata(imdep),ydata(jmdep)
181  REAL entree(imdep,jmdep)
182 c
183  INTEGER imar, jmar
184  REAL x(imar),y(jmar)
185  REAL sortie(imar,jmar)
186 c
187  INTEGER i, j, ii, jj
188  REAL a(400),b(400),c(200),d(200)
189  REAL number(400,200)
190  INTEGER indx(400,200), indy(400,200)
191  REAL dist(400,200), distsom(400,200)
192 c
193  IF (imar.GT.400 .OR. jmar.GT.200) THEN
194  print*, 'imar ou jmar trop grand', imar, jmar
195  CALL abort
196  ENDIF
197 c
198  IF (imdep.GT.400 .OR. jmdep.GT.200) THEN
199  print*, 'imdep ou jmdep trop grand', imdep, jmdep
200  CALL abort
201  ENDIF
202 c
203 c calculer les bords a et b de la nouvelle grille
204 c
205  a(1) = x(1) - (x(2)-x(1))/2.0
206  b(1) = (x(1)+x(2))/2.0
207  DO i = 2, imar-1
208  a(i) = b(i-1)
209  b(i) = (x(i)+x(i+1))/2.0
210  ENDDO
211  a(imar) = b(imar-1)
212  b(imar) = x(imar) + (x(imar)-x(imar-1))/2.0
213 
214 c
215 c calculer les bords c et d de la nouvelle grille
216 c
217  c(1) = y(1) - (y(2)-y(1))/2.0
218  d(1) = (y(1)+y(2))/2.0
219  DO j = 2, jmar-1
220  c(j) = d(j-1)
221  d(j) = (y(j)+y(j+1))/2.0
222  ENDDO
223  c(jmar) = d(jmar-1)
224  d(jmar) = y(jmar) + (y(jmar)-y(jmar-1))/2.0
225 
226 c
227 c trouver les indices (indx,indy) de la nouvelle grille sur laquelle
228 c un point de l'ancienne grille est tombe.
229 c
230 c
231 c ..... Modif P. Le Van ( 23/08/95 ) ....
232 
233  DO ii = 1, imar
234  DO jj = 1, jmar
235  DO i = 1, imdep
236  IF( ( xdata(i)-a(ii).GE.1.e-5.AND.xdata(i)-b(ii).LE.1.e-5 ).OR.
237  . ( xdata(i)-a(ii).LE.1.e-5.AND.xdata(i)-b(ii).GE.1.e-5 ) )
238  . THEN
239  DO j = 1, jmdep
240  IF( (ydata(j)-c(jj).GE.1.e-5.AND.ydata(j)-d(jj).LE.1.e-5 ).OR.
241  . ( ydata(j)-c(jj).LE.1.e-5.AND.ydata(j)-d(jj).GE.1.e-5 ) )
242  . THEN
243  indx(i,j) = ii
244  indy(i,j) = jj
245  ENDIF
246  ENDDO
247  ENDIF
248  ENDDO
249  ENDDO
250  ENDDO
251 c
252 c faire une verification
253 c
254 
255  DO i = 1, imdep
256  DO j = 1, jmdep
257  IF (indx(i,j).GT.imar .OR. indy(i,j).GT.jmar) THEN
258  print*, 'Probleme grave,i,j,indx,indy=',
259  . i,j,indx(i,j),indy(i,j)
260  CALL abort
261  ENDIF
262  ENDDO
263  ENDDO
264 
265 c
266 c calculer la distance des anciens points avec le nouveau point,
267 c on prend ensuite une sorte d'inverse pour ponderation.
268 c
269 
270  DO i = 1, imar
271  DO j = 1, jmar
272  number(i,j) = 0.0
273  distsom(i,j) = 0.0
274  ENDDO
275  ENDDO
276  DO i = 1, imdep
277  DO j = 1, jmdep
278  dist(i,j) = sqrt( (xdata(i)-x(indx(i,j)))**2
279  . +(ydata(j)-y(indy(i,j)))**2 )
280  distsom(indx(i,j),indy(i,j)) = distsom(indx(i,j),indy(i,j))
281  . + dist(i,j)
282  number(indx(i,j),indy(i,j)) = number(indx(i,j),indy(i,j)) +1.
283  ENDDO
284  ENDDO
285  DO i = 1, imdep
286  DO j = 1, jmdep
287  dist(i,j) = 1.0 - dist(i,j)/distsom(indx(i,j),indy(i,j))
288  ENDDO
289  ENDDO
290 
291  DO i = 1, imar
292  DO j = 1, jmar
293  number(i,j) = 0.0
294  sortie(i,j) = 0.0
295  ENDDO
296  ENDDO
297  DO i = 1, imdep
298  DO j = 1, jmdep
299  sortie(indx(i,j),indy(i,j)) = sortie(indx(i,j),indy(i,j))
300  . + entree(i,j) * dist(i,j)
301  number(indx(i,j),indy(i,j)) = number(indx(i,j),indy(i,j))
302  . + dist(i,j)
303  ENDDO
304  ENDDO
305  DO i = 1, imar
306  DO j = 1, jmar
307  IF (number(i,j) .GT. 0.001) THEN
308  sortie(i,j) = sortie(i,j) / number(i,j)
309  ELSE
310  print*, 'probleme,i,j=', i,j
311  CALL abort
312  ENDIF
313  ENDDO
314  ENDDO
315 
316  RETURN
317  END
318 
319 
320 
321 
322  SUBROUTINE mask_c_o(imdep, jmdep, xdata, ydata, relief,
323  . imar, jmar, x, y, mask)
324 c=======================================================================
325 c z.x.li (le 1 avril 1994): A partir du champ de relief, on fabrique
326 c un champ indicateur (masque) terre/ocean
327 c terre:1; ocean:0
328 c
329 c Methode naive (voir grille_m)
330 C=======================================================================
331  IMPLICIT none
332 
333  INTEGER imdep, jmdep
334  REAL xdata(imdep),ydata(jmdep)
335  REAL relief(imdep,jmdep)
336 c
337  INTEGER imar, jmar
338  REAL x(imar),y(jmar)
339  REAL mask(imar,jmar)
340 c
341  INTEGER i, j, ii, jj
342  REAL a(2200),b(2200),c(1100),d(1100)
343  REAL num_tot(2200,1100), num_oce(2200,1100)
344 c
345  IF (imar.GT.2200 .OR. jmar.GT.1100) THEN
346  print*, 'imar ou jmar trop grand', imar, jmar
347  CALL abort
348  ENDIF
349 c
350 
351  a(1) = x(1) - (x(2)-x(1))/2.0
352  b(1) = (x(1)+x(2))/2.0
353  DO i = 2, imar-1
354  a(i) = b(i-1)
355  b(i) = (x(i)+x(i+1))/2.0
356  ENDDO
357  a(imar) = b(imar-1)
358  b(imar) = x(imar) + (x(imar)-x(imar-1))/2.0
359 
360  c(1) = y(1) - (y(2)-y(1))/2.0
361  d(1) = (y(1)+y(2))/2.0
362  DO j = 2, jmar-1
363  c(j) = d(j-1)
364  d(j) = (y(j)+y(j+1))/2.0
365  ENDDO
366  c(jmar) = d(jmar-1)
367  d(jmar) = y(jmar) + (y(jmar)-y(jmar-1))/2.0
368 
369  DO i = 1, imar
370  DO j = 1, jmar
371  num_oce(i,j) = 0.0
372  num_tot(i,j) = 0.0
373  ENDDO
374  ENDDO
375 
376 c
377 c ..... Modif P. Le Van ( 23/08/95 ) ....
378 
379  DO ii = 1, imar
380  DO jj = 1, jmar
381  DO i = 1, imdep
382  IF( ( xdata(i)-a(ii).GE.1.e-5.AND.xdata(i)-b(ii).LE.1.e-5 ).OR.
383  . ( xdata(i)-a(ii).LE.1.e-5.AND.xdata(i)-b(ii).GE.1.e-5 ) )
384  . THEN
385  DO j = 1, jmdep
386  IF( (ydata(j)-c(jj).GE.1.e-5.AND.ydata(j)-d(jj).LE.1.e-5 ).OR.
387  . ( ydata(j)-c(jj).LE.1.e-5.AND.ydata(j)-d(jj).GE.1.e-5 ) )
388  . THEN
389  num_tot(ii,jj) = num_tot(ii,jj) + 1.0
390  IF (.NOT. ( relief(i,j) - 0.9. ge. 1.e-5 ) )
391  . num_oce(ii,jj) = num_oce(ii,jj) + 1.0
392  ENDIF
393  ENDDO
394  ENDIF
395  ENDDO
396  ENDDO
397  ENDDO
398 c
399 c
400 c
401  DO i = 1, imar
402  DO j = 1, jmar
403  IF (num_tot(i,j) .GT. 0.001) THEN
404  IF ( num_oce(i,j)/num_tot(i,j) - 0.5 .GE. 1.e-5 ) THEN
405  mask(i,j) = 0.
406  ELSE
407  mask(i,j) = 1.
408  ENDIF
409  ELSE
410  print*, 'probleme,i,j=', i,j
411  CALL abort
412  ENDIF
413  ENDDO
414  ENDDO
415 
416  RETURN
417  END
418 c
419 c
420 
421 
422  SUBROUTINE rugosite(imdep, jmdep, xdata, ydata, entree,
423  . imar, jmar, x, y, sortie, mask)
424 c=======================================================================
425 c z.x.li (le 1 avril 1994): Transformer la longueur de rugosite d'une
426 c grille fine a une grille grossiere. Sur l'ocean, on impose une valeur
427 c fixe (0.001m).
428 c
429 c Methode naive (voir grille_m)
430 C=======================================================================
431  IMPLICIT none
432 
433  INTEGER imdep, jmdep
434  REAL xdata(imdep),ydata(jmdep)
435  REAL entree(imdep,jmdep)
436 c
437  INTEGER imar, jmar
438  REAL x(imar),y(jmar)
439  REAL sortie(imar,jmar), mask(imar,jmar)
440 c
441  INTEGER i, j, ii, jj
442  REAL a(400),b(400),c(400),d(400)
443  REAL num_tot(400,400)
444  REAL distans(400*400)
445  INTEGER i_proche, j_proche, ij_proche
446 #ifdef CRAY
447  INTEGER ismin
448 #else
449  REAL zzmin
450 #endif
451 c
452  IF (imar.GT.400 .OR. jmar.GT.400) THEN
453  print*, 'imar ou jmar trop grand', imar, jmar
454  CALL abort
455  ENDIF
456 c
457 
458  a(1) = x(1) - (x(2)-x(1))/2.0
459  b(1) = (x(1)+x(2))/2.0
460  DO i = 2, imar-1
461  a(i) = b(i-1)
462  b(i) = (x(i)+x(i+1))/2.0
463  ENDDO
464  a(imar) = b(imar-1)
465  b(imar) = x(imar) + (x(imar)-x(imar-1))/2.0
466 
467  c(1) = y(1) - (y(2)-y(1))/2.0
468  d(1) = (y(1)+y(2))/2.0
469  DO j = 2, jmar-1
470  c(j) = d(j-1)
471  d(j) = (y(j)+y(j+1))/2.0
472  ENDDO
473  c(jmar) = d(jmar-1)
474  d(jmar) = y(jmar) + (y(jmar)-y(jmar-1))/2.0
475 
476  DO i = 1, imar
477  DO j = 1, jmar
478  num_tot(i,j) = 0.0
479  sortie(i,j) = 0.0
480  ENDDO
481  ENDDO
482 
483 c
484 c
485 c ..... Modif P. Le Van ( 23/08/95 ) ....
486 
487  DO ii = 1, imar
488  DO jj = 1, jmar
489  DO i = 1, imdep
490  IF( ( xdata(i)-a(ii).GE.1.e-5.AND.xdata(i)-b(ii).LE.1.e-5 ).OR.
491  . ( xdata(i)-a(ii).LE.1.e-5.AND.xdata(i)-b(ii).GE.1.e-5 ) )
492  . THEN
493  DO j = 1, jmdep
494  IF( (ydata(j)-c(jj).GE.1.e-5.AND.ydata(j)-d(jj).LE.1.e-5 ).OR.
495  . ( ydata(j)-c(jj).LE.1.e-5.AND.ydata(j)-d(jj).GE.1.e-5 ) )
496  . THEN
497  sortie(ii,jj) = sortie(ii,jj) + log(entree(i,j))
498  num_tot(ii,jj) = num_tot(ii,jj) + 1.0
499  ENDIF
500  ENDDO
501  ENDIF
502  ENDDO
503  ENDDO
504  ENDDO
505 c
506 
507  DO i = 1, imar
508  DO j = 1, jmar
509  IF (nint(mask(i,j)).EQ.1) THEN
510  IF (num_tot(i,j) .GT. 0.0) THEN
511  sortie(i,j) = sortie(i,j) / num_tot(i,j)
512  sortie(i,j) = exp(sortie(i,j))
513  ELSE
514  print*, 'probleme,i,j=', i,j
515 ccc CALL ABORT
516  CALL dist_sphe(x(i),y(j),xdata,ydata,imdep,jmdep,distans)
517 #ifdef CRAY
518  ij_proche = ismin(imdep*jmdep,distans,1)
519 #else
520  ij_proche = 1
521  zzmin = distans(ij_proche)
522  DO ii = 2, imdep*jmdep
523  IF (distans(ii).LT.zzmin) THEN
524  zzmin = distans(ii)
525  ij_proche = ii
526  ENDIF
527  ENDDO
528 #endif
529  j_proche = (ij_proche-1)/imdep + 1
530  i_proche = ij_proche - (j_proche-1)*imdep
531  print*, "solution:", ij_proche, i_proche, j_proche
532  sortie(i,j) = entree(i_proche,j_proche)
533  ENDIF
534  ELSE
535  sortie(i,j) = 0.001
536  ENDIF
537  ENDDO
538  ENDDO
539 
540  RETURN
541  END
542 
543 
544 
545 
546 
547  SUBROUTINE sea_ice(imdep, jmdep, xdata, ydata, glace01,
548  . imar, jmar, x, y, frac_ice)
549 c=======================================================================
550 c z.x.li (le 1 avril 1994): Transformer un champ d'indicateur de la
551 c glace (1, sinon 0) d'une grille fine a un champ de fraction de glace
552 c (entre 0 et 1) dans une grille plus grossiere.
553 c
554 c Methode naive (voir grille_m)
555 C=======================================================================
556  IMPLICIT none
557 
558  INTEGER imdep, jmdep
559  REAL xdata(imdep),ydata(jmdep)
560  REAL glace01(imdep,jmdep)
561 c
562  INTEGER imar, jmar
563  REAL x(imar),y(jmar)
564  REAL frac_ice(imar,jmar)
565 c
566  INTEGER i, j, ii, jj
567  REAL a(400),b(400),c(400),d(400)
568  REAL num_tot(400,400), num_ice(400,400)
569  REAL distans(400*400)
570  INTEGER i_proche, j_proche, ij_proche
571 #ifdef CRAY
572  INTEGER ismin
573 #else
574  REAL zzmin
575 #endif
576 c
577  IF (imar.GT.400 .OR. jmar.GT.400) THEN
578  print*, 'imar ou jmar trop grand', imar, jmar
579  CALL abort
580  ENDIF
581 c
582 
583  a(1) = x(1) - (x(2)-x(1))/2.0
584  b(1) = (x(1)+x(2))/2.0
585  DO i = 2, imar-1
586  a(i) = b(i-1)
587  b(i) = (x(i)+x(i+1))/2.0
588  ENDDO
589  a(imar) = b(imar-1)
590  b(imar) = x(imar) + (x(imar)-x(imar-1))/2.0
591 
592  c(1) = y(1) - (y(2)-y(1))/2.0
593  d(1) = (y(1)+y(2))/2.0
594  DO j = 2, jmar-1
595  c(j) = d(j-1)
596  d(j) = (y(j)+y(j+1))/2.0
597  ENDDO
598  c(jmar) = d(jmar-1)
599  d(jmar) = y(jmar) + (y(jmar)-y(jmar-1))/2.0
600 
601  DO i = 1, imar
602  DO j = 1, jmar
603  num_ice(i,j) = 0.0
604  num_tot(i,j) = 0.0
605  ENDDO
606  ENDDO
607 
608 c
609 c
610 c ..... Modif P. Le Van ( 23/08/95 ) ....
611 
612  DO ii = 1, imar
613  DO jj = 1, jmar
614  DO i = 1, imdep
615  IF( ( xdata(i)-a(ii).GE.1.e-5.AND.xdata(i)-b(ii).LE.1.e-5 ).OR.
616  . ( xdata(i)-a(ii).LE.1.e-5.AND.xdata(i)-b(ii).GE.1.e-5 ) )
617  . THEN
618  DO j = 1, jmdep
619  IF( (ydata(j)-c(jj).GE.1.e-5.AND.ydata(j)-d(jj).LE.1.e-5 ).OR.
620  . ( ydata(j)-c(jj).LE.1.e-5.AND.ydata(j)-d(jj).GE.1.e-5 ) )
621  . THEN
622  num_tot(ii,jj) = num_tot(ii,jj) + 1.0
623  IF (nint(glace01(i,j)).EQ.1 )
624  . num_ice(ii,jj) = num_ice(ii,jj) + 1.0
625  ENDIF
626  ENDDO
627  ENDIF
628  ENDDO
629  ENDDO
630  ENDDO
631 c
632 c
633 
634  DO i = 1, imar
635  DO j = 1, jmar
636  IF (num_tot(i,j) .GT. 0.001) THEN
637  IF (num_ice(i,j).GT.0.001) THEN
638  frac_ice(i,j) = num_ice(i,j) / num_tot(i,j)
639  ELSE
640  frac_ice(i,j) = 0.0
641  ENDIF
642  ELSE
643  print*, 'probleme,i,j=', i,j
644 ccc CALL ABORT
645  CALL dist_sphe(x(i),y(j),xdata,ydata,imdep,jmdep,distans)
646 #ifdef CRAY
647  ij_proche = ismin(imdep*jmdep,distans,1)
648 #else
649  ij_proche = 1
650  zzmin = distans(ij_proche)
651  DO ii = 2, imdep*jmdep
652  IF (distans(ii).LT.zzmin) THEN
653  zzmin = distans(ii)
654  ij_proche = ii
655  ENDIF
656  ENDDO
657 #endif
658  j_proche = (ij_proche-1)/imdep + 1
659  i_proche = ij_proche - (j_proche-1)*imdep
660  print*, "solution:", ij_proche, i_proche, j_proche
661  IF (nint(glace01(i_proche,j_proche)).EQ.1 ) THEN
662  frac_ice(i,j) = 1.0
663  ELSE
664  frac_ice(i,j) = 0.0
665  ENDIF
666  ENDIF
667  ENDDO
668  ENDDO
669 
670  RETURN
671  END
672 
673 
674 
675  SUBROUTINE rugsoro(imrel, jmrel, xrel, yrel, relief,
676  . immod, jmmod, xmod, ymod, rugs)
677 c=======================================================================
678 c Calculer la longueur de rugosite liee au relief en utilisant
679 c l'ecart-type dans une maille de 1x1
680 C=======================================================================
681  IMPLICIT none
682 c
683 #ifdef CRAY
684  INTEGER ismin
685 #else
686  REAL zzmin
687 #endif
688 c
689  REAL amin, amax
690 c
691  INTEGER imrel, jmrel
692  REAL xrel(imrel),yrel(jmrel)
693  REAL relief(imrel,jmrel)
694 c
695  INTEGER immod, jmmod
696  REAL xmod(immod),ymod(jmmod)
697  REAL rugs(immod,jmmod)
698 c
699  INTEGER imtmp, jmtmp
700  parameter(imtmp=360,jmtmp=180)
701  REAL xtmp(imtmp), ytmp(jmtmp)
702  REAL(KIND=8) cham1tmp(imtmp,jmtmp), cham2tmp(imtmp,jmtmp)
703  REAL zzzz
704 c
705  INTEGER i, j, ii, jj
706  REAL a(2200),b(2200),c(1100),d(1100)
707  REAL number(2200,1100)
708 c
709  REAL distans(400*400)
710  INTEGER i_proche, j_proche, ij_proche
711 c
712  IF (immod.GT.2200 .OR. jmmod.GT.1100) THEN
713  print*, 'immod ou jmmod trop grand', immod, jmmod
714  CALL abort
715  ENDIF
716 c
717 c Calculs intermediares:
718 c
719  xtmp(1) = -180.0 + 360.0/REAL(imtmp) / 2.0
720  DO i = 2, imtmp
721  xtmp(i) = xtmp(i-1) + 360.0/REAL(imtmp)
722  ENDDO
723  DO i = 1, imtmp
724  xtmp(i) = xtmp(i) /180.0 * 4.0*atan(1.0)
725  ENDDO
726  ytmp(1) = -90.0 + 180.0/REAL(jmtmp) / 2.0
727  DO j = 2, jmtmp
728  ytmp(j) = ytmp(j-1) + 180.0/REAL(jmtmp)
729  ENDDO
730  DO j = 1, jmtmp
731  ytmp(j) = ytmp(j) /180.0 * 4.0*atan(1.0)
732  ENDDO
733 c
734  a(1) = xtmp(1) - (xtmp(2)-xtmp(1))/2.0
735  b(1) = (xtmp(1)+xtmp(2))/2.0
736  DO i = 2, imtmp-1
737  a(i) = b(i-1)
738  b(i) = (xtmp(i)+xtmp(i+1))/2.0
739  ENDDO
740  a(imtmp) = b(imtmp-1)
741  b(imtmp) = xtmp(imtmp) + (xtmp(imtmp)-xtmp(imtmp-1))/2.0
742 
743  c(1) = ytmp(1) - (ytmp(2)-ytmp(1))/2.0
744  d(1) = (ytmp(1)+ytmp(2))/2.0
745  DO j = 2, jmtmp-1
746  c(j) = d(j-1)
747  d(j) = (ytmp(j)+ytmp(j+1))/2.0
748  ENDDO
749  c(jmtmp) = d(jmtmp-1)
750  d(jmtmp) = ytmp(jmtmp) + (ytmp(jmtmp)-ytmp(jmtmp-1))/2.0
751 
752  DO i = 1, imtmp
753  DO j = 1, jmtmp
754  number(i,j) = 0.0
755  cham1tmp(i,j) = 0.0
756  cham2tmp(i,j) = 0.0
757  ENDDO
758  ENDDO
759 c
760 c
761 c
762 c ..... Modif P. Le Van ( 23/08/95 ) ....
763 
764  DO ii = 1, imtmp
765  DO jj = 1, jmtmp
766  DO i = 1, imrel
767  IF( ( xrel(i)-a(ii).GE.1.e-5.AND.xrel(i)-b(ii).LE.1.e-5 ).OR.
768  . ( xrel(i)-a(ii).LE.1.e-5.AND.xrel(i)-b(ii).GE.1.e-5 ) )
769  . THEN
770  DO j = 1, jmrel
771  IF( (yrel(j)-c(jj).GE.1.e-5.AND.yrel(j)-d(jj).LE.1.e-5 ).OR.
772  . ( yrel(j)-c(jj).LE.1.e-5.AND.yrel(j)-d(jj).GE.1.e-5 ) )
773  . THEN
774  number(ii,jj) = number(ii,jj) + 1.0
775  cham1tmp(ii,jj) = cham1tmp(ii,jj) + relief(i,j)
776  cham2tmp(ii,jj) = cham2tmp(ii,jj)
777  . + relief(i,j)*relief(i,j)
778  ENDIF
779  ENDDO
780  ENDIF
781  ENDDO
782  ENDDO
783  ENDDO
784 c
785 c
786  DO i = 1, imtmp
787  DO j = 1, jmtmp
788  IF (number(i,j) .GT. 0.001) THEN
789  cham1tmp(i,j) = cham1tmp(i,j) / number(i,j)
790  cham2tmp(i,j) = cham2tmp(i,j) / number(i,j)
791  zzzz=cham2tmp(i,j)-cham1tmp(i,j)**2
792  if (zzzz .lt. 0.0) then
793  if (zzzz .gt. -7.5) then
794  zzzz = 0.0
795  print*,'Pb rugsoro, -7.5 < zzzz < 0, => zzz = 0.0'
796  else
797  stop 'Pb rugsoro, zzzz <-7.5'
798  endif
799  endif
800  cham2tmp(i,j) = sqrt(zzzz)
801  ELSE
802  print*, 'probleme,i,j=', i,j
803  CALL abort
804  ENDIF
805  ENDDO
806  ENDDO
807 c
808  amin = cham2tmp(1,1)
809  amax = cham2tmp(1,1)
810  DO j = 1, jmtmp
811  DO i = 1, imtmp
812  IF (cham2tmp(i,j).GT.amax) amax = cham2tmp(i,j)
813  IF (cham2tmp(i,j).LT.amin) amin = cham2tmp(i,j)
814  ENDDO
815  ENDDO
816  print*, 'Ecart-type 1x1:', amin, amax
817 c
818 c
819 c
820  a(1) = xmod(1) - (xmod(2)-xmod(1))/2.0
821  b(1) = (xmod(1)+xmod(2))/2.0
822  DO i = 2, immod-1
823  a(i) = b(i-1)
824  b(i) = (xmod(i)+xmod(i+1))/2.0
825  ENDDO
826  a(immod) = b(immod-1)
827  b(immod) = xmod(immod) + (xmod(immod)-xmod(immod-1))/2.0
828 
829  c(1) = ymod(1) - (ymod(2)-ymod(1))/2.0
830  d(1) = (ymod(1)+ymod(2))/2.0
831  DO j = 2, jmmod-1
832  c(j) = d(j-1)
833  d(j) = (ymod(j)+ymod(j+1))/2.0
834  ENDDO
835  c(jmmod) = d(jmmod-1)
836  d(jmmod) = ymod(jmmod) + (ymod(jmmod)-ymod(jmmod-1))/2.0
837 c
838  DO i = 1, immod
839  DO j = 1, jmmod
840  number(i,j) = 0.0
841  rugs(i,j) = 0.0
842  ENDDO
843  ENDDO
844 c
845 c
846 c
847 c ..... Modif P. Le Van ( 23/08/95 ) ....
848 
849  DO ii = 1, immod
850  DO jj = 1, jmmod
851  DO i = 1, imtmp
852  IF( ( xtmp(i)-a(ii).GE.1.e-5.AND.xtmp(i)-b(ii).LE.1.e-5 ).OR.
853  . ( xtmp(i)-a(ii).LE.1.e-5.AND.xtmp(i)-b(ii).GE.1.e-5 ) )
854  . THEN
855  DO j = 1, jmtmp
856  IF( (ytmp(j)-c(jj).GE.1.e-5.AND.ytmp(j)-d(jj).LE.1.e-5 ).OR.
857  . ( ytmp(j)-c(jj).LE.1.e-5.AND.ytmp(j)-d(jj).GE.1.e-5 ) )
858  . THEN
859  number(ii,jj) = number(ii,jj) + 1.0
860  rugs(ii,jj) = rugs(ii,jj)
861  . + log(max(0.001_8,cham2tmp(i,j)))
862  ENDIF
863  ENDDO
864  ENDIF
865  ENDDO
866  ENDDO
867  ENDDO
868 c
869 c
870  DO i = 1, immod
871  DO j = 1, jmmod
872  IF (number(i,j) .GT. 0.001) THEN
873  rugs(i,j) = rugs(i,j) / number(i,j)
874  rugs(i,j) = exp(rugs(i,j))
875  ELSE
876  print*, 'probleme,i,j=', i,j
877 ccc CALL ABORT
878  CALL dist_sphe(xmod(i),ymod(j),xtmp,ytmp,imtmp,jmtmp,distans)
879 #ifdef CRAY
880  ij_proche = ismin(imtmp*jmtmp,distans,1)
881 #else
882  ij_proche = 1
883  zzmin = distans(ij_proche)
884  DO ii = 2, imtmp*jmtmp
885  IF (distans(ii).LT.zzmin) THEN
886  zzmin = distans(ii)
887  ij_proche = ii
888  ENDIF
889  ENDDO
890 #endif
891  j_proche = (ij_proche-1)/imtmp + 1
892  i_proche = ij_proche - (j_proche-1)*imtmp
893  print*, "solution:", ij_proche, i_proche, j_proche
894  rugs(i,j) = log(max(0.001_8,cham2tmp(i_proche,j_proche)))
895  ENDIF
896  ENDDO
897  ENDDO
898 c
899  amin = rugs(1,1)
900  amax = rugs(1,1)
901  DO j = 1, jmmod
902  DO i = 1, immod
903  IF (rugs(i,j).GT.amax) amax = rugs(i,j)
904  IF (rugs(i,j).LT.amin) amin = rugs(i,j)
905  ENDDO
906  ENDDO
907  print*, 'Ecart-type du modele:', amin, amax
908 c
909  DO j = 1, jmmod
910  DO i = 1, immod
911  rugs(i,j) = rugs(i,j) / amax * 20.0
912  ENDDO
913  ENDDO
914 c
915  amin = rugs(1,1)
916  amax = rugs(1,1)
917  DO j = 1, jmmod
918  DO i = 1, immod
919  IF (rugs(i,j).GT.amax) amax = rugs(i,j)
920  IF (rugs(i,j).LT.amin) amin = rugs(i,j)
921  ENDDO
922  ENDDO
923  print*, 'Longueur de rugosite du modele:', amin, amax
924 c
925  RETURN
926  END
927 c
928  SUBROUTINE dist_sphe(rf_lon,rf_lat,rlon,rlat,im,jm,distance)
929 c
930 c Auteur: Laurent Li (le 30 decembre 1996)
931 c
932 c Ce programme calcule la distance minimale (selon le grand cercle)
933 c entre deux points sur la terre
934 c
935 c Input:
936  INTEGER im, jm ! dimensions
937  REAL rf_lon ! longitude du point de reference (degres)
938  REAL rf_lat ! latitude du point de reference (degres)
939  REAL rlon(im), rlat(jm) ! longitude et latitude des points
940 c
941 c Output:
942  REAL distance(im,jm) ! distances en metre
943 c
944  REAL rlon1, rlat1
945  REAL rlon2, rlat2
946  REAL dist
947  REAL pa, pb, p, pi
948 c
949  REAL radius
950  parameter(radius=6371229.)
951 c
952  pi = 4.0 * atan(1.0)
953 c
954  DO 9999 j = 1, jm
955  DO 9999 i = 1, im
956 c
957  rlon1=rf_lon
958  rlat1=rf_lat
959  rlon2=rlon(i)
960  rlat2=rlat(j)
961  pa = pi/2.0 - rlat1*pi/180.0 ! dist. entre pole n et point a
962  pb = pi/2.0 - rlat2*pi/180.0 ! dist. entre pole n et point b
963  p = (rlon1-rlon2)*pi/180.0 ! angle entre a et b (leurs meridiens)
964 c
965  dist = acos( cos(pa)*cos(pb) + sin(pa)*sin(pb)*cos(p))
966  dist = radius * dist
967  distance(i,j) = dist
968 c
969  9999 CONTINUE
970 c
971  END