Context Navigation

grid_atob.F @ 5434

Last change on this file since 5434 was 774, checked in by Laurent Fairhead, 18 years ago
Suite du merge entre la version et la HEAD: quelques modifications de Yann sur le LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 25.8 KB

Line
1	!
2	! $Header$
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*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/FLOAT(imtmp) / 2.0
720	DO i = 2, imtmp
721	xtmp(i) = xtmp(i-1) + 360.0/FLOAT(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/FLOAT(jmtmp) / 2.0
727	DO j = 2, jmtmp
728	ytmp(j) = ytmp(j-1) + 180.0/FLOAT(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,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,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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format