1 | PROGRAM xvik |
---|
2 | |
---|
3 | USE filtreg_mod, ONLY: inifilr |
---|
4 | USE comconst_mod, ONLY: dtvr,g,r,pi |
---|
5 | |
---|
6 | IMPLICIT NONE |
---|
7 | c======================================================================= |
---|
8 | c |
---|
9 | c Pression au site Viking |
---|
10 | c |
---|
11 | c======================================================================= |
---|
12 | c----------------------------------------------------------------------- |
---|
13 | c declarations: |
---|
14 | c ------------- |
---|
15 | |
---|
16 | |
---|
17 | #include "dimensions.h" |
---|
18 | #include "paramet.h" |
---|
19 | #include "comdissip.h" |
---|
20 | #include "comgeom2.h" |
---|
21 | !#include "control.h" |
---|
22 | #include "netcdf.inc" |
---|
23 | |
---|
24 | |
---|
25 | INTEGER itau,nbpas,nbpasmx |
---|
26 | PARAMETER(nbpasmx=1000000) |
---|
27 | REAL temps(nbpasmx) |
---|
28 | INTEGER unitlec |
---|
29 | INTEGER i,j,l,jj |
---|
30 | REAL constR |
---|
31 | |
---|
32 | c Declarations NCDF: |
---|
33 | c ----------------- |
---|
34 | CHARACTER*100 varname |
---|
35 | INTEGER ierr,nid,nvarid,dimid |
---|
36 | LOGICAL nc |
---|
37 | INTEGER start_ps(3),start_temp(4),start_co2ice(3) |
---|
38 | INTEGER count_ps(3),count_temp(4),count_co2ice(3) |
---|
39 | |
---|
40 | c declarations pour les points viking: |
---|
41 | c ------------------------------------ |
---|
42 | INTEGER ivik(2),jvik(2),ifile(2),iv |
---|
43 | REAL lonvik(2),latvik(2),phivik(2),phisim(2) |
---|
44 | REAL unanj |
---|
45 | |
---|
46 | c variables meteo: |
---|
47 | c ---------------- |
---|
48 | REAL vnat(iip1,jjm,llm),unat(iip1,jjp1,llm) |
---|
49 | REAL t(iip1,jjp1,llm),ps(iip1,jjp1),pstot, phis(iip1,jjp1) |
---|
50 | REAL co2ice(iip1,jjp1), captotN,captotS |
---|
51 | real t7(iip1,jjp1) ! temperature in 7th atmospheric layer |
---|
52 | |
---|
53 | REAL zp1,zp2,zp2_sm,zu,zv,zw(0:1,0:1,2),zalpha,zbeta |
---|
54 | |
---|
55 | LOGICAL firstcal,lcal,latcal,lvent,day_ls |
---|
56 | INTEGER*4 day0 |
---|
57 | |
---|
58 | REAL ziceco2(iip1,jjp1) |
---|
59 | REAL day,zt,sollong,sol,dayw |
---|
60 | REAL airtot1,gh |
---|
61 | |
---|
62 | INTEGER ii,iyear,kyear |
---|
63 | |
---|
64 | CHARACTER*2 chr2 |
---|
65 | |
---|
66 | |
---|
67 | c declarations de l'interface avec mywrite: |
---|
68 | c ----------------------------------------- |
---|
69 | |
---|
70 | CHARACTER file*80 |
---|
71 | CHARACTER pathchmp*80,pathsor*80,nomfich*80 |
---|
72 | |
---|
73 | c externe: |
---|
74 | c -------- |
---|
75 | |
---|
76 | EXTERNAL iniconst,inigeom,covcont,mywrite |
---|
77 | EXTERNAL exner,pbar |
---|
78 | EXTERNAL solarlong,coordij,moy2 |
---|
79 | EXTERNAL SSUM |
---|
80 | REAL SSUM |
---|
81 | |
---|
82 | cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
---|
83 | |
---|
84 | c----------------------------------------------------------------------- |
---|
85 | c initialisations: |
---|
86 | c ---------------- |
---|
87 | |
---|
88 | unanj=667.9 |
---|
89 | print*,'WARNING!!!',unanj,'Jours/an' |
---|
90 | nc=.true. |
---|
91 | lcal=.true. |
---|
92 | latcal=.true. |
---|
93 | lvent=.false. |
---|
94 | day_ls=.true. |
---|
95 | |
---|
96 | c lecture du fichier xvik.def |
---|
97 | |
---|
98 | phivik(1)=-3627 |
---|
99 | phivik(2)=-4505 |
---|
100 | |
---|
101 | |
---|
102 | |
---|
103 | OPEN(99,file='xvik.def',form='formatted') |
---|
104 | |
---|
105 | READ(99,*) |
---|
106 | READ(99,*,iostat=ierr) phivik |
---|
107 | IF(ierr.NE.0) GOTO 105 |
---|
108 | |
---|
109 | READ(99,*,END=105) |
---|
110 | READ(99,'(a)',END=105) pathchmp |
---|
111 | READ(99,*,END=105) |
---|
112 | READ(99,'(a)',END=105) pathsor |
---|
113 | READ(99,*,END=105) |
---|
114 | c READ(99,'(l1)',END=105) day_ls |
---|
115 | READ(99,'(l1)',END=105) |
---|
116 | READ(99,'(l1)',END=105) lcal |
---|
117 | READ(99,'(l1)',END=105) |
---|
118 | READ(99,'(l1)',END=105) lvent |
---|
119 | READ(99,'(l1)',END=105) |
---|
120 | READ(99,'(l1)',END=105) latcal |
---|
121 | |
---|
122 | 105 CONTINUE |
---|
123 | CLOSE(99) |
---|
124 | write (*,*)'>>>>>>>>>>>>>>>>', phivik,g |
---|
125 | DO iv=1,2 |
---|
126 | phivik(iv)=phivik(iv)*3.73 |
---|
127 | END DO |
---|
128 | |
---|
129 | write(*,*) ' pathchmp:',trim(pathchmp) |
---|
130 | write(*,*) ' pathsor:',trim(pathsor) |
---|
131 | |
---|
132 | c----------------------------------------------------------------------- |
---|
133 | c----------------------------------------------------------------------- |
---|
134 | c ouverture des fichiers xgraph: |
---|
135 | c ------------------------------ |
---|
136 | |
---|
137 | ifile(1)=12 |
---|
138 | ifile(2)=13 |
---|
139 | kyear=-1 |
---|
140 | c OPEN(77,file='xlongday',form='formatted') |
---|
141 | |
---|
142 | unitlec=11 |
---|
143 | c |
---|
144 | PRINT*,'entrer le nom du fichier NC' |
---|
145 | READ(5,'(a)') nomfich |
---|
146 | |
---|
147 | PRINT *,'nomfich : ',nomfich |
---|
148 | |
---|
149 | |
---|
150 | c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
151 | c grande boucle sur les fichiers histoire: |
---|
152 | c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
153 | |
---|
154 | firstcal=.true. |
---|
155 | DO WHILE(len_trim(nomfich).GT.0.AND.len_trim(nomfich).LT.50) |
---|
156 | PRINT *,'>>> nomfich : ',trim(nomfich) |
---|
157 | |
---|
158 | c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
159 | |
---|
160 | file=pathchmp(1:len_trim(pathchmp))//'/'// |
---|
161 | s nomfich(1:len_trim(nomfich)) |
---|
162 | PRINT*,'file.nc: ', file(1:len_trim(file))//'.nc' |
---|
163 | PRINT*,'timestep ',dtvr |
---|
164 | |
---|
165 | IF(nc) THEN |
---|
166 | ierr= NF_OPEN(file(1:len_trim(file))//'.nc',NF_NOWRITE,nid) |
---|
167 | ELSE |
---|
168 | PRINT*,'Ouverture binaire ',file |
---|
169 | OPEN(unitlec,file=file,status='old',form='unformatted', |
---|
170 | . iostat=ierr) |
---|
171 | ENDIF |
---|
172 | |
---|
173 | c---------------------------------------------------------------------- |
---|
174 | c initialisation de la physique: |
---|
175 | c ------------------------------ |
---|
176 | |
---|
177 | CALL readhead_NC(file(1:len_trim(file))//'.nc',day0,phis,constR) |
---|
178 | |
---|
179 | WRITE (*,*) 'day0 = ' , day0 |
---|
180 | |
---|
181 | CALL iniconst |
---|
182 | CALL inigeom |
---|
183 | CALL inifilr |
---|
184 | |
---|
185 | |
---|
186 | c Lecture temps : |
---|
187 | |
---|
188 | ierr= NF_INQ_DIMID (nid,"Time",dimid) |
---|
189 | IF (ierr.NE.NF_NOERR) THEN |
---|
190 | PRINT*, 'xvik: Le champ <Time> est absent' |
---|
191 | CALL abort |
---|
192 | ENDIF |
---|
193 | |
---|
194 | ierr= NF_INQ_DIMLEN (nid,dimid,nbpas) |
---|
195 | |
---|
196 | ierr = NF_INQ_VARID (nid, "Time", nvarid) |
---|
197 | #ifdef NC_DOUBLE |
---|
198 | ierr = NF_GET_VAR_DOUBLE(nid, nvarid, temps) |
---|
199 | #else |
---|
200 | ierr = NF_GET_VAR_REAL(nid, nvarid, temps) |
---|
201 | #endif |
---|
202 | IF (ierr.NE.NF_NOERR) THEN |
---|
203 | PRINT*, 'xvik: Lecture echouee pour <Time>' |
---|
204 | CALL abort |
---|
205 | ENDIF |
---|
206 | |
---|
207 | PRINT*,'temps',(temps(itau),itau=1,10) |
---|
208 | |
---|
209 | c----------------------------------------------------------------------- |
---|
210 | c coordonnees des point Viking: |
---|
211 | c ----------------------------- |
---|
212 | |
---|
213 | latvik(1)=22.27*pi/180. |
---|
214 | lonvik(1)=-47.9*pi/180. |
---|
215 | latvik(2)=47.67*pi/180. |
---|
216 | lonvik(2)=(360.-225.71)*pi/180. |
---|
217 | |
---|
218 | c ponderations pour les 4 points autour de Viking |
---|
219 | DO iv=1,2 |
---|
220 | CALL coordij(lonvik(iv),latvik(iv),ivik(iv),jvik(iv)) |
---|
221 | IF(lonvik(iv).lt.rlonv(ivik(iv))) THEN |
---|
222 | ivik(iv)=ivik(iv)-1 |
---|
223 | ENDIF |
---|
224 | IF(latvik(iv).gt.rlatu(jvik(iv))) THEN |
---|
225 | jvik(iv)=jvik(iv)-1 |
---|
226 | ENDIF |
---|
227 | zalpha=(lonvik(iv)-rlonv(ivik(iv)))/ |
---|
228 | s (rlonv(ivik(iv)+1)-rlonv(ivik(iv))) |
---|
229 | zbeta=(latvik(iv)-rlatu(jvik(iv)))/ |
---|
230 | s (rlatu(jvik(iv)+1)-rlatu(jvik(iv))) |
---|
231 | zw(0,0,iv)=(1.-zalpha)*(1.-zbeta) |
---|
232 | zw(1,0,iv)=zalpha*(1.-zbeta) |
---|
233 | zw(0,1,iv)=(1.-zalpha)*zbeta |
---|
234 | zw(1,1,iv)=zalpha*zbeta |
---|
235 | ENDDO |
---|
236 | |
---|
237 | c altitude reelle et modele aux points Viking |
---|
238 | DO iv=1,2 |
---|
239 | phisim(iv)=0. |
---|
240 | DO jj=0,1 |
---|
241 | j=jvik(iv)+jj |
---|
242 | DO ii=0,1 |
---|
243 | i=ivik(iv)+ii |
---|
244 | phisim(iv)=phisim(iv)+zw(ii,jj,iv)*phis(i,j) |
---|
245 | ENDDO |
---|
246 | ENDDO |
---|
247 | ENDDO |
---|
248 | PRINT*,'relief aux points Viking pour les sorties:',phivik |
---|
249 | |
---|
250 | c---------------------------------------------------------------------- |
---|
251 | c lectures des etats: |
---|
252 | c ------------------- |
---|
253 | |
---|
254 | airtot1=1./(SSUM(ip1jmp1,aire,1)-SSUM(jjp1,aire,iip1)) |
---|
255 | |
---|
256 | c====================================================================== |
---|
257 | c debut de la boucle sur les etats dans un fichier histoire: |
---|
258 | c====================================================================== |
---|
259 | count_ps=(/iip1,jjp1,1/) |
---|
260 | count_co2ice=(/iip1,jjp1,1/) |
---|
261 | count_temp=(/iip1,jjp1,llm,1/) |
---|
262 | |
---|
263 | DO itau=1,nbpas |
---|
264 | |
---|
265 | start_ps=(/1,1,itau/) |
---|
266 | start_co2ice=(/1,1,itau/) |
---|
267 | start_temp=(/1,1,1,itau/) |
---|
268 | c lecture drs des champs: |
---|
269 | c ----------------------- |
---|
270 | c varname='u' |
---|
271 | c ierr=drsread (unitlec,varname,unat,itau) |
---|
272 | c PRINT*,'unat',unat(iip1/2,jjp1/2,llm/2) |
---|
273 | c varname='v' |
---|
274 | c ierr=drsread (unitlec,varname,vnat,itau) |
---|
275 | c PRINT*,'vnat',vnat(iip1/2,jjp1/2,llm/2) |
---|
276 | |
---|
277 | ccccccccc LECTURE Ps ccccccccccccccccccccccccccc |
---|
278 | ierr = NF_INQ_VARID (nid, "ps", nvarid) |
---|
279 | #ifdef NC_DOUBLE |
---|
280 | ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start_ps,count_ps, ps) |
---|
281 | #else |
---|
282 | ierr = NF_GET_VARA_REAL(nid, nvarid,start_ps,count_ps, ps) |
---|
283 | #endif |
---|
284 | IF (ierr.NE.NF_NOERR) THEN |
---|
285 | PRINT*, 'xvik: Lecture echouee pour <ps>' |
---|
286 | CALL abort |
---|
287 | ENDIF |
---|
288 | |
---|
289 | PRINT*,'ps',ps(iip1/2,jjp1/2) |
---|
290 | |
---|
291 | ccccccccc LECTURE Temperature ccccccccccccccccccccccccccc |
---|
292 | ierr = NF_INQ_VARID (nid, "temp", nvarid) |
---|
293 | #ifdef NC_DOUBLE |
---|
294 | ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start_temp,count_temp, t) |
---|
295 | #else |
---|
296 | ierr = NF_GET_VARA_REAL(nid,nvarid,start_temp,count_temp, t) |
---|
297 | #endif |
---|
298 | IF (ierr.NE.NF_NOERR) THEN |
---|
299 | PRINT*, 'xvik: Lecture echouee pour <temp>' |
---|
300 | ! Ehouarn: proceed anyways |
---|
301 | ! CALL abort |
---|
302 | write(*,*)'--> Setting temperature to zero !!!' |
---|
303 | t(1:iip1,1:jjp1,1:llm)=0.0 |
---|
304 | write(*,*)'--> looking for temp7 (temp in 7th layer)' |
---|
305 | ierr=NF_INQ_VARID(nid,"temp7", nvarid) |
---|
306 | if (ierr.eq.NF_NOERR) then |
---|
307 | write(*,*) " OK, found temp7 variable" |
---|
308 | #ifdef NC_DOUBLE |
---|
309 | ierr=NF_GET_VARA_DOUBLE(nid,nvarid,start_ps,count_ps,t7) |
---|
310 | #else |
---|
311 | ierr=NF_GET_VARA_REAL(nid,nvarid,start_ps,count_ps,t7) |
---|
312 | #endif |
---|
313 | if (ierr.ne.NF_NOERR) then |
---|
314 | write(*,*)'xvik: failed loading temp7 !' |
---|
315 | stop |
---|
316 | endif |
---|
317 | else ! no 'temp7' variable |
---|
318 | write(*,*)' No temp7 variable either !' |
---|
319 | write(*,*)' Will have to to without ...' |
---|
320 | t7(1:iip1,1:jjp1)=0.0 |
---|
321 | endif |
---|
322 | ELSE ! t() was successfully loaded, copy 7th layer to t7() |
---|
323 | t7(1:iip1,1:jjp1)=t(1:iip1,1:jjp1,7) |
---|
324 | ENDIF |
---|
325 | |
---|
326 | c PRINT*,'t',t(iip1/2,jjp1/2,llm/2) |
---|
327 | |
---|
328 | ccccccccc LECTURE co2ice ccccccccccccccccccccccccccc |
---|
329 | ierr = NF_INQ_VARID (nid, "co2ice", nvarid) |
---|
330 | #ifdef NC_DOUBLE |
---|
331 | ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start_co2ice, |
---|
332 | & count_co2ice, co2ice) |
---|
333 | #else |
---|
334 | ierr = NF_GET_VARA_REAL(nid, nvarid,start_co2ice, |
---|
335 | & count_co2ice, co2ice) |
---|
336 | #endif |
---|
337 | IF (ierr.NE.NF_NOERR) THEN |
---|
338 | PRINT*, 'xvik: Lecture echouee pour <co2ice>' |
---|
339 | CALL abort |
---|
340 | ENDIF |
---|
341 | |
---|
342 | |
---|
343 | c Gestion du temps |
---|
344 | c ---------------- |
---|
345 | day=temps(itau) |
---|
346 | PRINT*,'day ',day |
---|
347 | CALL solarlong(day+day0,sollong) |
---|
348 | sol=day+day0+461. |
---|
349 | iyear=sol/unanj |
---|
350 | WRITE (*,*) 'iyear',iyear |
---|
351 | sol=sol-iyear*unanj |
---|
352 | c |
---|
353 | c Ouverture / fermeture des fichiers |
---|
354 | c ---------------------------------- |
---|
355 | IF (iyear.NE.kyear) THEN |
---|
356 | WRITE(chr2(1:1),'(i1)') iyear+1 |
---|
357 | WRITE (*,*) 'iyear bis',iyear |
---|
358 | WRITE (*,*) 'chr2' |
---|
359 | WRITE (*,*) chr2 |
---|
360 | IF(iyear.GE.9) WRITE(chr2,'(i2)') iyear+1 |
---|
361 | kyear=iyear |
---|
362 | DO ii=1,2 |
---|
363 | CLOSE(10+ifile(ii)) |
---|
364 | CLOSE(2+ifile(ii)) |
---|
365 | CLOSE(4+ifile(ii)) |
---|
366 | CLOSE(6+ifile(ii)) |
---|
367 | CLOSE(8+ifile(ii)) |
---|
368 | CLOSE(16+ifile(ii)) |
---|
369 | CLOSE(12+ifile(ii)) |
---|
370 | CLOSE(14+ifile(ii)) |
---|
371 | CLOSE(97) |
---|
372 | CLOSE(98) |
---|
373 | ENDDO |
---|
374 | CLOSE(5+ifile(1)) |
---|
375 | OPEN(ifile(1)+10,file='xpsol1'//chr2,form='formatted') |
---|
376 | OPEN(ifile(2)+10,file='xpsol2'//chr2,form='formatted') |
---|
377 | c OPEN(ifile(1)+8,file='xbpsol1'//chr2,form='formatted') |
---|
378 | c OPEN(ifile(2)+8,file='xbpsol2'//chr2,form='formatted') |
---|
379 | c OPEN(ifile(1)+2,file='xlps1'//chr2,form='formatted') |
---|
380 | c OPEN(ifile(2)+2,file='xlps2'//chr2,form='formatted') |
---|
381 | IF(lcal) THEN |
---|
382 | c OPEN(ifile(2)+4,file='xpressud'//chr2,form='formatted') |
---|
383 | c OPEN(ifile(1)+4,file='xpresnord'//chr2,form='formatted') |
---|
384 | c OPEN(ifile(1)+6,file='xpm2'//chr2,form='formatted') |
---|
385 | ENDIF |
---|
386 | IF(latcal) THEN |
---|
387 | c OPEN(ifile(2)+14,file='xlats'//chr2,form='formatted') |
---|
388 | c OPEN(ifile(1)+14,file='xlatn'//chr2,form='formatted') |
---|
389 | ENDIF |
---|
390 | IF(lvent) THEN |
---|
391 | c OPEN(ifile(1)+16,file='xu1'//chr2,form='formatted') |
---|
392 | c OPEN(ifile(2)+16,file='xu2'//chr2,form='formatted') |
---|
393 | c OPEN(ifile(1)+12,file='xv1'//chr2,form='formatted') |
---|
394 | c OPEN(ifile(2)+12,file='xv2'//chr2,form='formatted') |
---|
395 | ENDIF |
---|
396 | OPEN(97,file='xprestot'//chr2,form='formatted') |
---|
397 | c OPEN(98,file='xlat37_'//chr2,form='formatted') |
---|
398 | WRITE(98,'(f5.1,16f7.1)') 0.,(rlonv(i)*180./pi,i=1,iim,4) |
---|
399 | ENDIF |
---|
400 | |
---|
401 | |
---|
402 | sollong=sollong*180./pi |
---|
403 | IF(day_ls) THEN |
---|
404 | dayw=sol |
---|
405 | write(*,*) 'dayw', dayw |
---|
406 | ELSE |
---|
407 | dayw=sollong |
---|
408 | ENDIF |
---|
409 | |
---|
410 | c Calcul de la moyenne planetaire |
---|
411 | c ------------------------------- |
---|
412 | pstot=0. |
---|
413 | captotS=0. |
---|
414 | captotN=0. |
---|
415 | DO j=1,jjp1 |
---|
416 | DO i=1,iim |
---|
417 | pstot=pstot+aire(i,j)*ps(i,j) |
---|
418 | ENDDO |
---|
419 | ENDDO |
---|
420 | |
---|
421 | DO j=1,jjp1/2 |
---|
422 | DO i=1,iim |
---|
423 | captotN = captotN +aire(i,j)*co2ice(i,j) |
---|
424 | ENDDO |
---|
425 | ENDDO |
---|
426 | DO j=jjp1/2+1, jjp1 |
---|
427 | DO i=1,iim |
---|
428 | captotS = captotS +aire(i,j)*co2ice(i,j) |
---|
429 | ENDDO |
---|
430 | ENDDO |
---|
431 | WRITE(97,'(4e16.6)') dayw,pstot*airtot1 |
---|
432 | & , captotN*g*airtot1, captotS*g*airtot1 |
---|
433 | |
---|
434 | IF(.NOT.firstcal) THEN |
---|
435 | WRITE(98,'(f5.1,16f7.3)') |
---|
436 | s dayw,(ps(i,37),i=1,iim,4) |
---|
437 | |
---|
438 | c boucle sur les sites vikings: |
---|
439 | c ---------------------------- |
---|
440 | |
---|
441 | DO iv=1,2 |
---|
442 | |
---|
443 | c interpolation de la temperature dans la 7eme couche, de la pression |
---|
444 | c de surface et des vents aux points viking. |
---|
445 | |
---|
446 | zp1=0. |
---|
447 | zp2=0. |
---|
448 | zp2_sm=0. |
---|
449 | zt=0. |
---|
450 | zu=0. |
---|
451 | zv=0. |
---|
452 | DO jj=0,1 |
---|
453 | j=jvik(iv)+jj |
---|
454 | DO ii=0,1 |
---|
455 | i=ivik(iv)+ii |
---|
456 | ! zt=zt+zw(ii,jj,iv)*t(i,j,7) |
---|
457 | zt=zt+zw(ii,jj,iv)*t7(i,j) |
---|
458 | ! zp1=zp1+zw(ii,jj,iv)*ps(i,j) |
---|
459 | zp1=zp1+zw(ii,jj,iv)*log(ps(i,j)) ! interpolate in log(P) |
---|
460 | WRITE (*,*) 'ps autour iv',ps(i,j),iv |
---|
461 | zu=zu+zw(ii,jj,iv)*unat(i,j,1)/cu(i,j) |
---|
462 | zv=zv+zw(ii,jj,iv)*vnat(i,j,1)/cv(i,j) |
---|
463 | ENDDO |
---|
464 | ENDDO |
---|
465 | zp1=exp(zp1) ! because of the bilinear interpolation in log(P) |
---|
466 | |
---|
467 | c pression au sol extrapolee a partir de la temp. 7eme couche |
---|
468 | WRITE (*,*) 'constR ',constR |
---|
469 | WRITE (*,*) 'zt ',zt |
---|
470 | gh=constR*zt |
---|
471 | if (gh.eq.0) then ! if we don't have temperature values |
---|
472 | ! assume a scale height of 10km |
---|
473 | zp2=zp1*exp(-(phivik(iv)-phisim(iv))/(3.73*1.e4)) |
---|
474 | else |
---|
475 | zp2=zp1*exp(-(phivik(iv)-phisim(iv))/gh) |
---|
476 | endif |
---|
477 | WRITE (*,*) 'iv,pstot,zp2, zp1, phivik(iv),phisim(iv),gh' |
---|
478 | WRITE (*,*) iv,pstot*airtot1,zp2,zp1,phivik(iv),phisim(iv),gh |
---|
479 | ! WRITE(ifile(iv)+10,'(2e15.5)') dayw,zp1 |
---|
480 | WRITE(ifile(iv)+10,'(3e15.5)') dayw,zp2,zp1 |
---|
481 | |
---|
482 | c sorties eventuelles de vent |
---|
483 | IF(lvent) THEN |
---|
484 | WRITE(ifile(iv)+16,'(2e15.5)') |
---|
485 | s dayw,zu |
---|
486 | WRITE(ifile(iv)+12,'(2e15.5)') |
---|
487 | s dayw,zv |
---|
488 | ENDIF |
---|
489 | ENDDO |
---|
490 | c IF (lcal) THEN |
---|
491 | c WRITE(ifile(1)+4,'(2e15.6)') dayw,airtot1*g*.01* |
---|
492 | c s (SSUM(ip1jmp1/2,ziceco2,1)-SSUM(jjp1/2,ziceco2,iip1)) |
---|
493 | c WRITE(ifile(2)+4,'(2e15.6)') dayw,airtot1*g*.01* |
---|
494 | c s (SSUM(iip1*jjm/2,ziceco2(1,jjm/2+2),1)- |
---|
495 | c s SSUM(jjm/2,ziceco2(1,jjm/2+2),iip1)) |
---|
496 | c ENDIF |
---|
497 | c IF(latcal) THEN |
---|
498 | c CALL icelat(iim,jjm,ziceco2,rlatv,zicelat) |
---|
499 | c WRITE(ifile(1)+14,'(2e15.6)') dayw,zicelat(1)*180./pi |
---|
500 | c WRITE(ifile(2)+14,'(2e15.6)') dayw,zicelat(2)*180./pi |
---|
501 | c ENDIF |
---|
502 | ENDIF |
---|
503 | firstcal=.false. |
---|
504 | |
---|
505 | c====================================================================== |
---|
506 | c Fin de la boucle sur les etats du fichier histoire: |
---|
507 | c====================================================================== |
---|
508 | ENDDO |
---|
509 | |
---|
510 | ierr= NF_CLOSE(nid) |
---|
511 | |
---|
512 | PRINT*,'Fin du fichier',nomfich |
---|
513 | print*,'Entrer un nouveau fichier ou return pour finir' |
---|
514 | READ(5,'(a)',err=9999) nomfich |
---|
515 | |
---|
516 | c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
517 | c Fin de la boucle sur les fichiers histoire: |
---|
518 | c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
519 | |
---|
520 | ENDDO |
---|
521 | |
---|
522 | PRINT*,'relief du point V1',.001*phis(ivik(1),jvik(1))/g |
---|
523 | PRINT*,'relief du point V2',.001*phis(ivik(2),jvik(2))/g |
---|
524 | DO iv=1,2 |
---|
525 | PRINT*,'Viking',iv,' i=',ivik(iv),'j =',jvik(iv) |
---|
526 | WRITE(6,7777) |
---|
527 | s (rlonv(i)*180./pi,i=ivik(iv)-1,ivik(iv)+2) |
---|
528 | print* |
---|
529 | DO j=jvik(iv)-1,jvik(iv)+2 |
---|
530 | WRITE(6,'(f8.1,10x,5f7.1)') |
---|
531 | s rlatu(j)*180./pi,(phis(i,j)/(g*1000.),i=ivik(iv)-1,ivik(iv)+2) |
---|
532 | ENDDO |
---|
533 | print* |
---|
534 | print*,'zw' |
---|
535 | write(6,'(2(2f10.4/))') ((zw(ii,jj,iv),ii=0,1),jj=0,1) |
---|
536 | print*,'altitude interpolee (km) ',phisim(iv)/1000./g |
---|
537 | ENDDO |
---|
538 | PRINT*,'R=',r |
---|
539 | 9999 PRINT*,'Fin ' |
---|
540 | |
---|
541 | 7777 FORMAT ('latitude/longitude',4f7.1) |
---|
542 | END |
---|