source: trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/xvik.F @ 1564

Last change on this file since 1564 was 1422, checked in by milmd, 10 years ago

In GENERIC, MARS and COMMON models replace some include files by modules (usefull for decoupling physics with dynamics).

File size: 16.8 KB
Line 
1      PROGRAM xvik
2
3      USE filtreg_mod, ONLY: inifilr
4      USE comconst_mod, ONLY: dtvr,g,r,pi
5
6      IMPLICIT NONE
7c=======================================================================
8c
9c  Pression au site Viking
10c
11c=======================================================================
12c-----------------------------------------------------------------------
13c   declarations:
14c   -------------
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
32c   Declarations NCDF:
33c   -----------------
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
40c   declarations pour les points viking:
41c   ------------------------------------
42      INTEGER ivik(2),jvik(2),ifile(2),iv
43      REAL lonvik(2),latvik(2),phivik(2),phisim(2)
44      REAL unanj
45
46c   variables meteo:
47c   ----------------
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       
67c   declarations de l'interface avec mywrite:
68c   -----------------------------------------
69
70      CHARACTER file*80
71      CHARACTER pathchmp*80,pathsor*80,nomfich*80
72
73c   externe:
74c   --------
75
76      EXTERNAL iniconst,inigeom,covcont,mywrite
77      EXTERNAL exner,pbar
78      EXTERNAL solarlong,coordij,moy2
79      EXTERNAL SSUM
80      REAL SSUM
81
82cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
83
84c-----------------------------------------------------------------------
85c   initialisations:
86c   ----------------
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
96c 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)
114c     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     
132c-----------------------------------------------------------------------
133c-----------------------------------------------------------------------
134c   ouverture des fichiers xgraph:
135c   ------------------------------
136
137      ifile(1)=12
138      ifile(2)=13
139      kyear=-1
140c      OPEN(77,file='xlongday',form='formatted')
141
142      unitlec=11
143c
144      PRINT*,'entrer le nom du fichier NC'
145      READ(5,'(a)') nomfich
146
147      PRINT *,'nomfich : ',nomfich
148
149
150c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
151c   grande boucle sur les fichiers histoire:
152c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
153
154      firstcal=.true.
155      DO WHILE(len_trim(nomfich).GT.0.AND.len_trim(nomfich).LT.50)
156      PRINT *,'>>>  nomfich : ',trim(nomfich)
157
158c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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
173c----------------------------------------------------------------------
174c   initialisation de la physique:
175c   ------------------------------
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
186c   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             
209c-----------------------------------------------------------------------
210c   coordonnees des point Viking:
211c   -----------------------------
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
218c   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
237c   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
250c----------------------------------------------------------------------
251c   lectures des etats:
252c   -------------------
253
254       airtot1=1./(SSUM(ip1jmp1,aire,1)-SSUM(jjp1,aire,iip1))
255
256c======================================================================
257c   debut de la boucle sur les etats dans un fichier histoire:
258c======================================================================
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/)
268c   lecture drs des champs:
269c   -----------------------
270c         varname='u'
271c         ierr=drsread (unitlec,varname,unat,itau)
272c         PRINT*,'unat',unat(iip1/2,jjp1/2,llm/2)
273c         varname='v'
274c         ierr=drsread (unitlec,varname,vnat,itau)
275c         PRINT*,'vnat',vnat(iip1/2,jjp1/2,llm/2)
276
277ccccccccc  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
291ccccccccc  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
326c          PRINT*,'t',t(iip1/2,jjp1/2,llm/2)
327
328ccccccccc  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
343c Gestion du temps
344c ----------------
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
352c
353c Ouverture / fermeture des fichiers
354c ----------------------------------
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')
377c            OPEN(ifile(1)+8,file='xbpsol1'//chr2,form='formatted')
378c            OPEN(ifile(2)+8,file='xbpsol2'//chr2,form='formatted')
379c            OPEN(ifile(1)+2,file='xlps1'//chr2,form='formatted')
380c            OPEN(ifile(2)+2,file='xlps2'//chr2,form='formatted')
381             IF(lcal) THEN
382c               OPEN(ifile(2)+4,file='xpressud'//chr2,form='formatted')
383c               OPEN(ifile(1)+4,file='xpresnord'//chr2,form='formatted')
384c               OPEN(ifile(1)+6,file='xpm2'//chr2,form='formatted')
385             ENDIF
386                         IF(latcal) THEN
387c               OPEN(ifile(2)+14,file='xlats'//chr2,form='formatted')
388c               OPEN(ifile(1)+14,file='xlatn'//chr2,form='formatted')
389                         ENDIF
390             IF(lvent) THEN
391c               OPEN(ifile(1)+16,file='xu1'//chr2,form='formatted')
392c               OPEN(ifile(2)+16,file='xu2'//chr2,form='formatted')
393c               OPEN(ifile(1)+12,file='xv1'//chr2,form='formatted')
394c               OPEN(ifile(2)+12,file='xv2'//chr2,form='formatted')
395             ENDIF
396             OPEN(97,file='xprestot'//chr2,form='formatted')
397c            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
410c Calcul de la moyenne planetaire
411c -------------------------------
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
438c boucle sur les sites vikings:
439c ----------------------------
440
441          DO iv=1,2
442
443c interpolation de la temperature dans la 7eme couche, de la pression
444c 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 
467c               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           
482c   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
490c         IF (lcal) THEN
491c            WRITE(ifile(1)+4,'(2e15.6)') dayw,airtot1*g*.01*
492c    s       (SSUM(ip1jmp1/2,ziceco2,1)-SSUM(jjp1/2,ziceco2,iip1))
493c            WRITE(ifile(2)+4,'(2e15.6)') dayw,airtot1*g*.01*
494c    s       (SSUM(iip1*jjm/2,ziceco2(1,jjm/2+2),1)-
495c    s       SSUM(jjm/2,ziceco2(1,jjm/2+2),iip1))
496c         ENDIF
497c            IF(latcal) THEN
498c               CALL icelat(iim,jjm,ziceco2,rlatv,zicelat)
499c               WRITE(ifile(1)+14,'(2e15.6)') dayw,zicelat(1)*180./pi
500c               WRITE(ifile(2)+14,'(2e15.6)') dayw,zicelat(2)*180./pi
501c            ENDIF
502         ENDIF
503         firstcal=.false.
504
505c======================================================================
506c   Fin de la boucle sur les etats du fichier histoire:
507c======================================================================
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
516c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
517c   Fin de la boucle sur les fichiers histoire:
518c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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
5417777  FORMAT ('latitude/longitude',4f7.1)
542      END
Note: See TracBrowser for help on using the repository browser.