/[lmdze]/trunk/dyn3d/grille_m.f
ViewVC logotype

Annotation of /trunk/dyn3d/grille_m.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 15 - (hide annotations)
Fri Aug 1 15:24:12 2008 UTC (15 years, 9 months ago) by guez
Original Path: trunk/libf/dyn3d/grid_atob.f90
File size: 27809 byte(s) 
-- Minor modification of input/output:

Added variable "Sigma_O3_Royer" to "histday.nc". "ecrit_day" is not
modified in "physiq". Removed variables "pyu1", "pyv1", "ftsol1",
"ftsol2", "ftsol3", "ftsol4", "psrf1", "psrf2", "psrf3", "psrf4"
"mfu", "mfd", "en_u", "en_d", "de_d", "de_u", "coefh" from
"histrac.nc".

Variable "raz_date" of module "conf_gcm_m" has logical type instead of
integer type.

-- Should not change any result at run time:

Modified calls to "IOIPSL_Lionel" procedures because the interfaces of
these procedures have been simplified.

Changed name of variable in module "start_init_orog_m": "masque" to
"mask".

Created a module containing procedure "phyredem".

Removed arguments "punjours", "pdayref" and "ptimestep" of procedure
"iniphysiq".

Renamed procedure "gr_phy_write" to "gr_phy_write_2d". Created
procedure "gr_phy_write_3d".

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