Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

grid_atob.F @ 937

Last change on this file since 937 was 38, checked in by emillour, 14 years ago
Ajout du modè Martien (mon LMDZ.MARS.BETA, du 28/01/2011) dans le rértoire mars, pour pouvoir suivre plus facilement les modifs. EM
File size: 25.8 KB

Line
1	SUBROUTINE grille_m(imdep, jmdep, xdata, ydata, entree,
2	. imar, jmar, x, y, sortie)
3	c=======================================================================
4	c z.x.li (le 1 avril 1994) (voir aussi A. Harzallah et L. Fairhead)
5	c
6	c Methode naive pour transformer un champ d'une grille fine a une
7	c grille grossiere. Je considere que les nouveaux points occupent
8	c une zone adjacente qui comprend un ou plusieurs anciens points
9	c
10	c Aucune ponderation est consideree (voir grille_p)
11	c
12	c (c)
13	c ----d-----
14	c \| . . . .\|
15	c \| \|
16	c (b)a . * . .b(a)
17	c \| \|
18	c \| . . . .\|
19	c ----c-----
20	c (d)
21	C=======================================================================
22	c INPUT:
23	c imdep, jmdep: dimensions X et Y pour depart
24	c xdata, ydata: coordonnees X et Y pour depart
25	c entree: champ d'entree a transformer
26	c OUTPUT:
27	c imar, jmar: dimensions X et Y d'arrivee
28	c x, y: coordonnees X et Y d'arrivee
29	c sortie: champ de sortie deja transforme
30	C=======================================================================
31	IMPLICIT none
32
33	INTEGER imdep, jmdep
34	REAL xdata(imdep),ydata(jmdep)
35	REAL entree(imdep,jmdep)
36	c
37	INTEGER imar, jmar
38	REAL x(imar),y(jmar)
39	REAL sortie(imar,jmar)
40	c
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
51	c
52	IF (imar.GT.2200 .OR. jmar.GT.1100) THEN
53	PRINT*, 'imar ou jmar trop grand', imar, jmar
54	CALL ABORT
55	ENDIF
56	c
57	c Calculer les limites des zones des nouveaux points
58	c
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
84	c
85	c Determiner la zone sur laquelle chaque ancien point se trouve
86	c
87	c
88	c ..... 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
108	c
109	c Si aucun ancien point tombe sur une zone, c'est un probleme
110	c
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
118	ccc 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)
146	c=======================================================================
147	c z.x.li (le 1 avril 1994) (voir aussi A. Harzallah et L. Fairhead)
148	c
149	c Methode naive pour transformer un champ d'une grille fine a une
150	c grille grossiere. Je considere que les nouveaux points occupent
151	c une zone adjacente qui comprend un ou plusieurs anciens points
152	c
153	c Consideration de la distance des points (voir grille_m)
154	c
155	c (c)
156	c ----d-----
157	c \| . . . .\|
158	c \| \|
159	c (b)a . * . .b(a)
160	c \| \|
161	c \| . . . .\|
162	c ----c-----
163	c (d)
164	C=======================================================================
165	c INPUT:
166	c imdep, jmdep: dimensions X et Y pour depart
167	c xdata, ydata: coordonnees X et Y pour depart
168	c entree: champ d'entree a transformer
169	c OUTPUT:
170	c imar, jmar: dimensions X et Y d'arrivee
171	c x, y: coordonnees X et Y d'arrivee
172	c sortie: champ de sortie deja transforme
173	C=======================================================================
174	IMPLICIT none
175
176	INTEGER imdep, jmdep
177	REAL xdata(imdep),ydata(jmdep)
178	REAL entree(imdep,jmdep)
179	c
180	INTEGER imar, jmar
181	REAL x(imar),y(jmar)
182	REAL sortie(imar,jmar)
183	c
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)
189	c
190	IF (imar.GT.400 .OR. jmar.GT.200) THEN
191	PRINT*, 'imar ou jmar trop grand', imar, jmar
192	CALL ABORT
193	ENDIF
194	c
195	IF (imdep.GT.400 .OR. jmdep.GT.200) THEN
196	PRINT*, 'imdep ou jmdep trop grand', imdep, jmdep
197	CALL ABORT
198	ENDIF
199	c
200	c calculer les bords a et b de la nouvelle grille
201	c
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
211	c
212	c calculer les bords c et d de la nouvelle grille
213	c
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	c
224	c trouver les indices (indx,indy) de la nouvelle grille sur laquelle
225	c un point de l'ancienne grille est tombe.
226	c
227	c
228	c ..... 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
248	c
249	c faire une verification
250	c
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
262	c
263	c calculer la distance des anciens points avec le nouveau point,
264	c on prend ensuite une sorte d'inverse pour ponderation.
265	c
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)
321	c=======================================================================
322	c z.x.li (le 1 avril 1994): A partir du champ de relief, on fabrique
323	c un champ indicateur (masque) terre/ocean
324	c terre:1; ocean:0
325	c
326	c Methode naive (voir grille_m)
327	C=======================================================================
328	IMPLICIT none
329
330	INTEGER imdep, jmdep
331	REAL xdata(imdep),ydata(jmdep)
332	REAL relief(imdep,jmdep)
333	c
334	INTEGER imar, jmar
335	REAL x(imar),y(jmar)
336	REAL mask(imar,jmar)
337	c
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)
341	c
342	IF (imar.GT.2200 .OR. jmar.GT.1100) THEN
343	PRINT*, 'imar ou jmar trop grand', imar, jmar
344	CALL ABORT
345	ENDIF
346	c
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
373	c
374	c ..... 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
395	c
396	c
397	c
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
415	c
416	c
417
418
419	SUBROUTINE rugosite(imdep, jmdep, xdata, ydata, entree,
420	. imar, jmar, x, y, sortie, mask)
421	c=======================================================================
422	c z.x.li (le 1 avril 1994): Transformer la longueur de rugosite d'une
423	c grille fine a une grille grossiere. Sur l'ocean, on impose une valeur
424	c fixe (0.001m).
425	c
426	c Methode naive (voir grille_m)
427	C=======================================================================
428	IMPLICIT none
429
430	INTEGER imdep, jmdep
431	REAL xdata(imdep),ydata(jmdep)
432	REAL entree(imdep,jmdep)
433	c
434	INTEGER imar, jmar
435	REAL x(imar),y(jmar)
436	REAL sortie(imar,jmar), mask(imar,jmar)
437	c
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
448	c
449	IF (imar.GT.400 .OR. jmar.GT.400) THEN
450	PRINT*, 'imar ou jmar trop grand', imar, jmar
451	CALL ABORT
452	ENDIF
453	c
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
480	c
481	c
482	c ..... 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
502	c
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
512	ccc 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)
546	c=======================================================================
547	c z.x.li (le 1 avril 1994): Transformer un champ d'indicateur de la
548	c glace (1, sinon 0) d'une grille fine a un champ de fraction de glace
549	c (entre 0 et 1) dans une grille plus grossiere.
550	c
551	c Methode naive (voir grille_m)
552	C=======================================================================
553	IMPLICIT none
554
555	INTEGER imdep, jmdep
556	REAL xdata(imdep),ydata(jmdep)
557	REAL glace01(imdep,jmdep)
558	c
559	INTEGER imar, jmar
560	REAL x(imar),y(jmar)
561	REAL frac_ice(imar,jmar)
562	c
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
573	c
574	IF (imar.GT.400 .OR. jmar.GT.400) THEN
575	PRINT*, 'imar ou jmar trop grand', imar, jmar
576	CALL ABORT
577	ENDIF
578	c
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
605	c
606	c
607	c ..... 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
628	c
629	c
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
641	ccc 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)
674	c=======================================================================
675	c Calculer la longueur de rugosite liee au relief en utilisant
676	c l'ecart-type dans une maille de 1x1
677	C=======================================================================
678	IMPLICIT none
679	c
680	#ifdef CRAY
681	INTEGER ISMIN
682	#else
683	REAL zzmin
684	#endif
685	c
686	REAL amin, AMAX
687	c
688	INTEGER imrel, jmrel
689	REAL xrel(imrel),yrel(jmrel)
690	REAL relief(imrel,jmrel)
691	c
692	INTEGER immod, jmmod
693	REAL xmod(immod),ymod(jmmod)
694	REAL rugs(immod,jmmod)
695	c
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
701	c
702	INTEGER i, j, ii, jj
703	REAL a(2200),b(2200),c(1100),d(1100)
704	REAL number(2200,1100)
705	c
706	REAL distans(400*400)
707	INTEGER i_proche, j_proche, ij_proche
708	c
709	IF (immod.GT.2200 .OR. jmmod.GT.1100) THEN
710	PRINT*, 'immod ou jmmod trop grand', immod, jmmod
711	CALL ABORT
712	ENDIF
713	c
714	c Calculs intermediares:
715	c
716	xtmp(1) = -180.0 + 360.0/REAL(imtmp) / 2.0
717	DO i = 2, imtmp
718	xtmp(i) = xtmp(i-1) + 360.0/REAL(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/REAL(jmtmp) / 2.0
724	DO j = 2, jmtmp
725	ytmp(j) = ytmp(j-1) + 180.0/REAL(jmtmp)
726	ENDDO
727	DO j = 1, jmtmp
728	ytmp(j) = ytmp(j) /180.0 * 4.0*ATAN(1.0)
729	ENDDO
730	c
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
756	c
757	c
758	c
759	c ..... 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
781	c
782	c
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
804	c
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
814	c
815	c
816	c
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
834	c
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
841	c
842	c
843	c
844	c ..... 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,real(cham2tmp(i,j))))
859	ENDIF
860	ENDDO
861	ENDIF
862	ENDDO
863	ENDDO
864	ENDDO
865	c
866	c
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
874	ccc 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,real(cham2tmp(i_proche,j_proche))))
892	ENDIF
893	ENDDO
894	ENDDO
895	c
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
905	c
906	DO j = 1, jmmod
907	DO i = 1, immod
908	rugs(i,j) = rugs(i,j) / AMAX * 20.0
909	ENDDO
910	ENDDO
911	c
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
921	c
922	RETURN
923	END
924	c
925	SUBROUTINE dist_sphe(rf_lon,rf_lat,rlon,rlat,im,jm,distance)
926	c
927	c Auteur: Laurent Li (le 30 decembre 1996)
928	c
929	c Ce programme calcule la distance minimale (selon le grand cercle)
930	c entre deux points sur la terre
931	c
932	c 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
937	c
938	c Output:
939	REAL distance(im,jm) ! distances en metre
940	c
941	REAL rlon1, rlat1
942	REAL rlon2, rlat2
943	REAL dist
944	REAL pa, pb, p, pi
945	c
946	REAL radius
947	PARAMETER (radius=6371229.)
948	c
949	pi = 4.0 * ATAN(1.0)
950	c
951	DO 9999 j = 1, jm
952	DO 9999 i = 1, im
953	c
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)
961	c
962	dist = ACOS( COS(pa)COS(pb) + SIN(pa)SIN(pb)*COS(p))
963	dist = radius * dist
964	distance(i,j) = dist
965	c
966	9999 CONTINUE
967	c
968	END

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

Download in other formats:

Original Format