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

Last change on this file since 1242 was 1130, checked in by emillour, 11 years ago

Mars GCM:
Series of changes to enable running in parallel (using LMDZ.COMMON dynamics);
Current LMDZ.MARS can still notheless be compiled and run in serial mode
"as previously".
Summary of main changes:

  • Main programs (newstart, start2archive, xvik) that used to be in dyn3d have been moved to phymars.
  • dyn3d/control.h is now module control_mod.F90
  • rearanged input/outputs routines everywhere to handle serial/MPI cases. physdem.F => phyredem.F90 , phyetat0.F => phyetat0.F90 ; all read/write routines for startfi files are gathered in module iostart.F90
  • added parallelism related routines init_phys_lmdz.F90, comgeomphy.F90, dimphy.F90, iniphysiq.F90, mod_grid_phy_lmdz.F90, mod_phys_lmdz_mpi_data.F90, mod_phys_lmdz_mpi_transfert.F90, mod_phys_lmdz_omp_data.F90, mod_phys_lmdz_omp_transfert.F90, mod_phys_lmdz_para.F90, mod_phys_lmdz_transfert_para.F90 in phymars and mod_const_mpi.F90 in dyn3d (for compliance with parallel case)
  • created generic routines 'planetwide_maxval' and 'planetwide_minval', in module "planetwide_mod", that enable obtaining the min and max of a field over the whole planet.

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