source: trunk/mesoscale/LMDZ.MARS/libf_gcm/dyn3d/xvik.F @ 57

Last change on this file since 57 was 57, checked in by aslmd, 14 years ago

mineur LMD_MM_MARS: ajout du GCM ancienne physique, systeme maintenant complet sur SVN (ne manque que la base de donnees d'etats initiaux)

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