source: LMDZ.3.3/branches/rel-LF/libf/dyn3d/grid_atob.F

Last change on this file was 2, checked in by lmdz, 25 years ago

Initial revision

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