Context Navigation

grid_atob.F @ 1108

Last change on this file since 1108 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

Rev	Line
[630]	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