source: LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/grid_atob.F @ 1322

Last change on this file since 1322 was 1299, checked in by Laurent Fairhead, 15 years ago

Nettoyage general pour se rapprocher des normes et éviter des erreurs a la
compilation:

  • tous les FLOAT() sont remplacés par des REAL()
  • tous les STOP dans phylmd sont remplacés par des appels à abort_gcm
  • le common control défini dans le fichier control.h est remplacé par le module control_mod pour éviter des messages sur l'alignement des variables dans les déclarations
  • des $Header$ remplacés par des $Id$ pour svn

Quelques remplacements à faire ont pu m'échapper


General cleanup of the code to try and adhere to norms and to prevent some
compilation errors:

  • all FLOAT() instructions have been replaced by REAL() instructions
  • all STOP instructions in phylmd have been replaced by calls to abort_gcm
  • the common block control defined in the control.h file has been replaced by the control_mod to prevent compilation warnings on the alignement of declared variables
  • $Header$ replaced by $Id$ for svn

Some changes which should have been made might have escaped me

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 25.8 KB
Line 
1!
2! $Id: grid_atob.F 1299 2010-01-20 14:27:21Z fairhead $
3!
4      SUBROUTINE grille_m(imdep, jmdep, xdata, ydata, entree,
5     .                    imar, jmar, x, y, sortie)
6c=======================================================================
7c z.x.li (le 1 avril 1994) (voir aussi A. Harzallah et L. Fairhead)
8c
9c Methode naive pour transformer un champ d'une grille fine a une
10c grille grossiere. Je considere que les nouveaux points occupent
11c une zone adjacente qui comprend un ou plusieurs anciens points
12c
13c Aucune ponderation est consideree (voir grille_p)
14c
15c           (c)
16c        ----d-----
17c        | . . . .|
18c        |        |
19c     (b)a . * . .b(a)
20c        |        |
21c        | . . . .|
22c        ----c-----
23c           (d)
24C=======================================================================
25c INPUT:
26c        imdep, jmdep: dimensions X et Y pour depart
27c        xdata, ydata: coordonnees X et Y pour depart
28c        entree: champ d'entree a transformer
29c OUTPUT:
30c        imar, jmar: dimensions X et Y d'arrivee
31c        x, y: coordonnees X et Y d'arrivee
32c        sortie: champ de sortie deja transforme
33C=======================================================================
34      IMPLICIT none
35
36      INTEGER imdep, jmdep
37      REAL xdata(imdep),ydata(jmdep)
38      REAL entree(imdep,jmdep)
39c
40      INTEGER imar, jmar
41      REAL x(imar),y(jmar)
42      REAL sortie(imar,jmar)
43c
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
54c
55      IF (imar.GT.2200 .OR. jmar.GT.1100) THEN
56         PRINT*, 'imar ou jmar trop grand', imar, jmar
57         CALL ABORT
58      ENDIF
59c
60c Calculer les limites des zones des nouveaux points
61c
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
87c
88c Determiner la zone sur laquelle chaque ancien point se trouve
89c
90c
91c  .....  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
111c
112c Si aucun ancien point tombe sur une zone, c'est un probleme
113c
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
121ccc         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)
149c=======================================================================
150c z.x.li (le 1 avril 1994) (voir aussi A. Harzallah et L. Fairhead)
151c
152c Methode naive pour transformer un champ d'une grille fine a une
153c grille grossiere. Je considere que les nouveaux points occupent
154c une zone adjacente qui comprend un ou plusieurs anciens points
155c
156c Consideration de la distance des points (voir grille_m)
157c
158c           (c)
159c        ----d-----
160c        | . . . .|
161c        |        |
162c     (b)a . * . .b(a)
163c        |        |
164c        | . . . .|
165c        ----c-----
166c           (d)
167C=======================================================================
168c INPUT:
169c        imdep, jmdep: dimensions X et Y pour depart
170c        xdata, ydata: coordonnees X et Y pour depart
171c        entree: champ d'entree a transformer
172c OUTPUT:
173c        imar, jmar: dimensions X et Y d'arrivee
174c        x, y: coordonnees X et Y d'arrivee
175c        sortie: champ de sortie deja transforme
176C=======================================================================
177      IMPLICIT none
178
179      INTEGER imdep, jmdep
180      REAL xdata(imdep),ydata(jmdep)
181      REAL entree(imdep,jmdep)
182c
183      INTEGER imar, jmar
184      REAL x(imar),y(jmar)
185      REAL sortie(imar,jmar)
186c
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)
192c
193      IF (imar.GT.400 .OR. jmar.GT.200) THEN
194         PRINT*, 'imar ou jmar trop grand', imar, jmar
195         CALL ABORT
196      ENDIF
197c
198      IF (imdep.GT.400 .OR. jmdep.GT.200) THEN
199         PRINT*, 'imdep ou jmdep trop grand', imdep, jmdep
200         CALL ABORT
201      ENDIF
202c
203c calculer les bords a et b de la nouvelle grille
204c
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
214c
215c calculer les bords c et d de la nouvelle grille
216c
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
226c
227c trouver les indices (indx,indy) de la nouvelle grille sur laquelle
228c un point de l'ancienne grille est tombe.
229c
230c
231c  .....  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
251c
252c faire une verification
253c
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
265c
266c calculer la distance des anciens points avec le nouveau point,
267c on prend ensuite une sorte d'inverse pour ponderation.
268c
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)
324c=======================================================================
325c z.x.li (le 1 avril 1994): A partir du champ de relief, on fabrique
326c                           un champ indicateur (masque) terre/ocean
327c                           terre:1; ocean:0
328c
329c Methode naive (voir grille_m)
330C=======================================================================
331      IMPLICIT none
332
333      INTEGER imdep, jmdep
334      REAL xdata(imdep),ydata(jmdep)
335      REAL relief(imdep,jmdep)
336c
337      INTEGER imar, jmar
338      REAL x(imar),y(jmar)
339      REAL mask(imar,jmar)
340c
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)
344c
345      IF (imar.GT.2200 .OR. jmar.GT.1100) THEN
346         PRINT*, 'imar ou jmar trop grand', imar, jmar
347         CALL ABORT
348      ENDIF
349c
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
376c
377c  .....  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
398c
399c
400c
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
418c
419c
420
421
422      SUBROUTINE rugosite(imdep, jmdep, xdata, ydata, entree,
423     .                    imar, jmar, x, y, sortie, mask)
424c=======================================================================
425c z.x.li (le 1 avril 1994): Transformer la longueur de rugosite d'une
426c grille fine a une grille grossiere. Sur l'ocean, on impose une valeur
427c fixe (0.001m).
428c
429c Methode naive (voir grille_m)
430C=======================================================================
431      IMPLICIT none
432
433      INTEGER imdep, jmdep
434      REAL xdata(imdep),ydata(jmdep)
435      REAL entree(imdep,jmdep)
436c
437      INTEGER imar, jmar
438      REAL x(imar),y(jmar)
439      REAL sortie(imar,jmar), mask(imar,jmar)
440c
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
451c
452      IF (imar.GT.400 .OR. jmar.GT.400) THEN
453         PRINT*, 'imar ou jmar trop grand', imar, jmar
454         CALL ABORT
455      ENDIF
456c
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
483c
484c
485c  .....  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
505c
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
515ccc            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)
549c=======================================================================
550c z.x.li (le 1 avril 1994): Transformer un champ d'indicateur de la
551c glace (1, sinon 0) d'une grille fine a un champ de fraction de glace
552c (entre 0 et 1) dans une grille plus grossiere.
553c
554c Methode naive (voir grille_m)
555C=======================================================================
556      IMPLICIT none
557
558      INTEGER imdep, jmdep
559      REAL xdata(imdep),ydata(jmdep)
560      REAL glace01(imdep,jmdep)
561c
562      INTEGER imar, jmar
563      REAL x(imar),y(jmar)
564      REAL frac_ice(imar,jmar)
565c
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
576c
577      IF (imar.GT.400 .OR. jmar.GT.400) THEN
578         PRINT*, 'imar ou jmar trop grand', imar, jmar
579         CALL ABORT
580      ENDIF
581c
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
608c
609c
610c  .....  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
631c
632c
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
644ccc           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)
677c=======================================================================
678c Calculer la longueur de rugosite liee au relief en utilisant
679c l'ecart-type dans une maille de 1x1
680C=======================================================================
681      IMPLICIT none
682c
683#ifdef CRAY
684      INTEGER ISMIN
685#else
686      REAL zzmin
687#endif
688c
689      REAL amin, AMAX
690c
691      INTEGER imrel, jmrel
692      REAL xrel(imrel),yrel(jmrel)
693      REAL relief(imrel,jmrel)
694c
695      INTEGER immod, jmmod
696      REAL xmod(immod),ymod(jmmod)
697      REAL rugs(immod,jmmod)
698c
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
704c
705      INTEGER i, j, ii, jj
706      REAL a(2200),b(2200),c(1100),d(1100)
707      REAL number(2200,1100)
708c
709      REAL distans(400*400)
710      INTEGER i_proche, j_proche, ij_proche
711c
712      IF (immod.GT.2200 .OR. jmmod.GT.1100) THEN
713         PRINT*, 'immod ou jmmod trop grand', immod, jmmod
714         CALL ABORT
715      ENDIF
716c
717c Calculs intermediares:
718c
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
733c
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
759c
760c
761c
762c  .....  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
784c
785c
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
807c
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
817c
818c
819c
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
837c
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
844c
845c
846c
847c  .....  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
868c
869c
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
877ccc         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
898c
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
908c
909      DO j = 1, jmmod
910      DO i = 1, immod
911         rugs(i,j) = rugs(i,j) / AMAX * 20.0
912      ENDDO
913      ENDDO
914c
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
924c
925      RETURN
926      END
927c
928      SUBROUTINE dist_sphe(rf_lon,rf_lat,rlon,rlat,im,jm,distance)
929c
930c Auteur: Laurent Li (le 30 decembre 1996)
931c
932c Ce programme calcule la distance minimale (selon le grand cercle)
933c entre deux points sur la terre
934c
935c 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
940c
941c Output:
942      REAL distance(im,jm) ! distances en metre
943c
944      REAL rlon1, rlat1
945      REAL rlon2, rlat2
946      REAL dist
947      REAL pa, pb, p, pi
948c
949      REAL radius
950      PARAMETER (radius=6371229.)
951c
952      pi = 4.0 * ATAN(1.0)
953c
954      DO 9999 j = 1, jm
955      DO 9999 i = 1, im
956c
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)
964c
965      dist = ACOS( COS(pa)*COS(pb) + SIN(pa)*SIN(pb)*COS(p))
966      dist = radius * dist
967      distance(i,j) = dist
968c
969 9999 CONTINUE
970c
971      END
Note: See TracBrowser for help on using the repository browser.