/[lmdze]/trunk/libf/dyn3d/grid_atob.f90
ViewVC logotype

Contents of /trunk/libf/dyn3d/grid_atob.f90

Parent Directory Parent Directory | Revision Log Revision Log

Revision 15 - (show annotations)
Fri Aug 1 15:24:12 2008 UTC (15 years, 8 months ago) by guez
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 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 function grille_m(xdata, ydata, entree, x, y)
10
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 use numer_rec, only: assert_eq
41
42 REAL, intent(in):: xdata(:),ydata(:)
43 REAL, intent(in):: entree(:, :)
44 REAL, intent(in):: x(:), y(:)
45
46 real grille_m(size(x), size(y))
47
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 use numer_rec, only: assert_eq
417
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 use numer_rec, only: assert_eq
536
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