source: trunk/LMDZ.TITAN/libf/phytitan/newstart.F @ 1056

Last change on this file since 1056 was 1056, checked in by slebonnois, 11 years ago

SL: Titan runs ! see DOC/chantiers/commit_importants.log

File size: 37.1 KB
Line 
1C======================================================================
2      PROGRAM newstart
3c=======================================================================
4c
5c
6c   Auteur:   S. Lebonnois,
7c    a partir des newstart/start_archive/lect_start_archive martiens
8c
9c             Derniere modif : 02/09 (ecriture des q*)
10c                              01/12 (inclusion dans svn dyn3d)
11c
12c   Objet:  Modify the grid for the initial state (LMD GCM VENUS/TITAN)
13c   -----           (from file NetCDF start_archive.nc)
14c
15c
16c=======================================================================
17
18      use IOIPSL
19      USE filtreg_mod
20      USE startvar
21      USE control_mod
22      USE infotrac
23      use cpdet_mod, only: ini_cpdet,t2tpot
24
25      implicit none
26
27#include "dimensions.h"
28#include "paramet.h"
29#include "comconst.h"
30#include "comdissnew.h"
31#include "comvert.h"
32#include "comgeom2.h"
33#include "logic.h"
34#include "temps.h"
35#include "ener.h"
36#include "description.h"
37#include "serre.h"
38#include "dimsoil.h"
39#include "netcdf.inc"
40
41c-----------------------------------------------------------------------
42c   Declarations
43c-----------------------------------------------------------------------
44
45c Variables pour fichier "ini"
46c------------------------------------
47      INTEGER   imold,jmold,lmold,nqold,ip1jmp1old
48      INTEGER   length
49      parameter (length = 100)
50      real      tab_cntrl(2*length)
51      INTEGER isoil,iq,iqmax
52      CHARACTER*2   str2
53
54c Variable histoire
55c------------------
56      REAL vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm) ! vents covariants
57      REAL teta(iip1,jjp1,llm),ps(iip1,jjp1)
58      REAL phis(iip1,jjp1)                     ! geopotentiel au sol
59      REAL masse(ip1jmp1,llm)                ! masse de l'atmosphere
60      REAL, ALLOCATABLE, DIMENSION(:,:,:,:):: q! champs advectes
61      REAL tab_cntrl_dyn(length) ! tableau des parametres de start
62
63c variable physique
64c------------------
65      integer    ngridmx
66      parameter (ngridmx=(2+(jjm-1)*iim - 1/jjm))
67      REAL tab_cntrl_fi(length) ! tableau des parametres de startfi
68      real rlat(ngridmx),rlon(ngridmx)
69      REAL tsurf(ngridmx),tsoil(ngridmx,nsoilmx)
70      REAL albe(ngridmx),radsol(ngridmx),sollw(ngridmx)
71      real solsw(ngridmx),dlw(ngridmx)
72      REAL zmea(ngridmx), zstd(ngridmx)
73      REAL zsig(ngridmx), zgam(ngridmx), zthe(ngridmx)
74      REAL zpic(ngridmx), zval(ngridmx)
75      real t_fi(ngridmx,llm)
76
77c Variable nouvelle grille naturelle au point scalaire
78c------------------------------------------------------
79      real us(iip1,jjp1,llm),vs(iip1,jjp1,llm)
80      REAL p3d(iip1,jjp1,llm+1)            ! pression aux interfaces
81      REAL phisold_newgrid(iip1,jjp1)
82      REAL T(iip1,jjp1,llm)
83      real rlonS(iip1,jjp1),rlatS(iip1,jjp1)
84      real tsurfS(iip1,jjp1),tsoilS(iip1,jjp1,nsoilmx)
85      real albeS(ip1jmp1),radsolS(ip1jmp1),sollwS(ip1jmp1)
86      real solswS(ip1jmp1),dlwS(ip1jmp1)
87      real zmeaS(ip1jmp1),zstdS(ip1jmp1),zsigS(ip1jmp1)
88      real zgamS(ip1jmp1),ztheS(ip1jmp1),zpicS(ip1jmp1)
89      real zvalS(ip1jmp1)
90
91      real ptotal
92
93c Var intermediaires : vent naturel, mais pas coord scalaire
94c-----------------------------------------------------------
95      real vnat(iip1,jjm,llm),unat(iip1,jjp1,llm)
96
97      REAL pks(iip1,jjp1)                      ! exner (f pour filtre)
98      REAL pk(iip1,jjp1,llm)
99      REAL pkf(iip1,jjp1,llm)
100      REAL alpha(iip1,jjp1,llm),beta(iip1,jjp1,llm)
101
102
103c Variable de l'ancienne grille
104c---------------------------------------------------------
105
106      real, dimension(:),     allocatable :: rlonuold, rlatvold
107      real, dimension(:),     allocatable :: rlonvold, rlatuold
108      real, dimension(:),     allocatable :: nivsigsold,nivsigold
109      real, dimension(:),     allocatable :: apold,bpold
110      real, dimension(:),     allocatable :: presnivsold
111      real, dimension(:,:,:), allocatable :: uold,vold,told
112      real, dimension(:,:,:,:), allocatable :: qold
113      real, dimension(:,:,:), allocatable :: tsoilold
114      real, dimension(:,:),   allocatable :: psold,phisold
115      real, dimension(:,:),   allocatable :: tsurfold
116      real, dimension(:,:),   allocatable :: albeold,radsolold
117      real, dimension(:,:),   allocatable :: sollwold,solswold
118      real, dimension(:,:),   allocatable :: dlwold
119      real, dimension(:,:),   allocatable :: zmeaold,zstdold,zsigold
120      real, dimension(:,:),   allocatable :: zgamold,ztheold,zpicold
121      real, dimension(:,:),   allocatable :: zvalold
122
123      real ptotalold
124
125c Variable intermediaires iutilise pour l'extrapolation verticale
126c----------------------------------------------------------------
127      real, dimension(:,:,:), allocatable :: var,varp1
128
129c divers local
130c-----------------
131
132      integer ierr,nid,nvarid
133      INTEGER ij, l,i,j
134      character*80      fichnom     
135      integer, dimension(4) :: start,counter
136      REAL phisinverse(iip1,jjp1)  ! geopotentiel au sol avant inversion
137      logical topoflag,albedoflag
138      real    albedo
139     
140c=======================================================================
141c  INITIALISATIONS DIVERSES
142c=======================================================================
143
144c VENUS/TITAN
145
146        iflag_trac = 1
147c-----------------------------------------------------------------------
148c   Initialisation des traceurs
149c   ---------------------------
150c  Choix du nombre de traceurs et du schema pour l'advection
151c  dans fichier traceur.def, par default ou via INCA
152      call infotrac_init
153
154c Allocation de la tableau q : champs advectes   
155      allocate(q(iip1,jjp1,llm,nqtot))
156
157c-----------------------------------------------------------------------
158c   Ouverture du fichier a modifier (start_archive.nc)
159c-----------------------------------------------------------------------
160
161        write(*,*) 'Creation d un etat initial a partir de'
162        write(*,*) './start_archive.nc'
163        write(*,*)
164        fichnom = 'start_archive.nc'
165        ierr = NF_OPEN (fichnom, NF_NOWRITE,nid)
166        IF (ierr.NE.NF_NOERR) THEN
167          write(6,*)' Pb d''ouverture du fichier ',fichnom
168          write(6,*)' ierr = ', ierr
169          CALL ABORT
170        ENDIF
171 
172c-----------------------------------------------------------------------
173c Lecture du tableau des parametres du run (pour la dynamique)
174c-----------------------------------------------------------------------
175
176        write(*,*) 'lecture tab_cntrl START_ARCHIVE'
177c
178        ierr = NF_INQ_VARID (nid, "controle", nvarid)
179#ifdef NC_DOUBLE
180        ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tab_cntrl)
181#else
182        ierr = NF_GET_VAR_REAL(nid, nvarid, tab_cntrl)
183#endif
184c
185      write(*,*) 'Impression de tab_cntrl'
186      do i=1,200
187        write(*,*) i,tab_cntrl(i)
188      enddo
189     
190c-----------------------------------------------------------------------
191c               Initialisation des constantes
192c-----------------------------------------------------------------------
193
194      imold      = tab_cntrl(1)
195      jmold      = tab_cntrl(2)
196      lmold      = tab_cntrl(3)
197      day_ref    = tab_cntrl(4)
198      annee_ref  = tab_cntrl(5)
199      rad        = tab_cntrl(6)
200      omeg       = tab_cntrl(7)
201      g          = tab_cntrl(8)
202      cpp        = tab_cntrl(9)
203      kappa      = tab_cntrl(10)
204      daysec     = tab_cntrl(11)
205      dtvr       = tab_cntrl(12)
206      etot0      = tab_cntrl(13)
207      ptot0      = tab_cntrl(14)
208      ztot0      = tab_cntrl(15)
209      stot0      = tab_cntrl(16)
210      ang0       = tab_cntrl(17)
211      pa         = tab_cntrl(18)
212      preff      = tab_cntrl(19)
213c
214      clon       = tab_cntrl(20)
215      clat       = tab_cntrl(21)
216      grossismx  = tab_cntrl(22)
217      grossismy  = tab_cntrl(23)
218
219      IF ( tab_cntrl(24).EQ.1. )  THEN
220        fxyhypb  = . TRUE .
221        dzoomx   = tab_cntrl(25)
222        dzoomy   = tab_cntrl(26)
223        taux     = tab_cntrl(28)
224        tauy     = tab_cntrl(29)
225      ELSE
226        fxyhypb = . FALSE .
227        ysinus  = . FALSE .
228        IF( tab_cntrl(27).EQ.1. ) ysinus = . TRUE.
229      ENDIF
230
231      ptotalold  = tab_cntrl(2*length)
232
233      write(*,*) "Old dimensions:"
234      write(*,*) "longitude: ",imold
235      write(*,*) "latitude: ",jmold
236      write(*,*) "altitude: ",lmold
237      write(*,*)
238
239      ip1jmp1old = (imold+1-1/imold)*(jmold+1-1/jmold)
240     
241c dans run.def
242      CALL getin('planet_type',planet_type)
243      call ini_cpdet
244
245c=======================================================================
246c   CHANGEMENT DE CONSTANTES CONTENUES DANS tab_cntrl
247c=======================================================================
248c  changement de la resolution dans le fichier de controle
249      tab_cntrl(1)  = REAL(iim)
250      tab_cntrl(2)  = REAL(jjm)
251      tab_cntrl(3)  = REAL(llm)
252
253c--------------------------------
254c We use a specific run.def to read new parameters that need to be changed
255c--------------------------------
256     
257c Changement de cpp:
258      CALL getin('cpp',cpp)
259      kappa = (8.314511/43.44E-3)/cpp
260      tab_cntrl(9)  = cpp
261      tab_cntrl(10) = kappa
262
263c CHANGING THE ZOOM PARAMETERS TO CHANGE THE GRID
264      CALL getin('clon',clon)
265      CALL getin('clat',clat)
266      tab_cntrl(20) = clon
267      tab_cntrl(21) = clat
268      CALL getin('grossismx',grossismx)
269      CALL getin('grossismy',grossismy)
270      tab_cntrl(22) = grossismx
271      tab_cntrl(23) = grossismy
272      CALL getin('fxyhypb',fxyhypb)
273      IF ( fxyhypb )  THEN
274        CALL getin('dzoomx',dzoomx)
275        CALL getin('dzoomy',dzoomy)
276        tab_cntrl(25) = dzoomx
277        tab_cntrl(26) = dzoomy
278        CALL getin('taux',taux)
279        CALL getin('tauy',tauy)
280        tab_cntrl(28) = taux
281        tab_cntrl(29) = tauy
282      ELSE
283        CALL getin('ysinus',ysinus)
284        tab_cntrl(27) = ysinus
285      ENDIF
286
287c changes must be done BEFORE these lines...
288      DO l=1,length
289         tab_cntrl_dyn(l)= tab_cntrl(l)
290         tab_cntrl_fi(l) = tab_cntrl(l+length)
291      ENDDO
292
293c-----------------------------------------------------------------------
294c       Autres initialisations
295c-----------------------------------------------------------------------
296
297      CALL iniconst
298      CALL inigeom
299      call inifilr
300
301c-----------------------------------------------------------------------
302c               Allocations des anciennes variables
303c-----------------------------------------------------------------------
304
305      allocate(rlonuold(imold+1), rlatvold(jmold))
306      allocate(rlonvold(imold+1), rlatuold(jmold+1))
307      allocate(nivsigsold(lmold+1),nivsigold(lmold))
308      allocate(apold(lmold),bpold(lmold))
309      allocate(presnivsold(lmold))
310      allocate(uold(imold+1,jmold+1,lmold))
311      allocate(vold(imold+1,jmold+1,lmold))
312      allocate(told(imold+1,jmold+1,lmold))
313      allocate(qold(imold+1,jmold+1,lmold,nqtot))
314      allocate(psold(imold+1,jmold+1))
315      allocate(phisold(imold+1,jmold+1))
316      allocate(tsurfold(imold+1,jmold+1))
317      allocate(tsoilold(imold+1,jmold+1,nsoilmx))
318      allocate(albeold(imold+1,jmold+1),radsolold(imold+1,jmold+1))
319      allocate(sollwold(imold+1,jmold+1),solswold(imold+1,jmold+1))
320      allocate(dlwold(imold+1,jmold+1))
321      allocate(zmeaold(imold+1,jmold+1),zstdold(imold+1,jmold+1))
322      allocate(zsigold(imold+1,jmold+1),zgamold(imold+1,jmold+1))
323      allocate(ztheold(imold+1,jmold+1),zpicold(imold+1,jmold+1))
324      allocate(zvalold(imold+1,jmold+1))
325
326      allocate(var (imold+1,jmold+1,llm))
327      allocate(varp1 (imold+1,jmold+1,llm+1))
328
329      print*,"Initialisations OK"
330
331c-----------------------------------------------------------------------
332c Lecture des longitudes et latitudes
333c-----------------------------------------------------------------------
334c
335      ierr = NF_INQ_VARID (nid, "rlonv", nvarid)
336      IF (ierr .NE. NF_NOERR) THEN
337         PRINT*, "new_grid: Le champ <rlonv> est absent de start.nc"
338         CALL abort
339      ENDIF
340#ifdef NC_DOUBLE
341      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlonvold)
342#else
343      ierr = NF_GET_VAR_REAL(nid, nvarid, rlonvold)
344#endif
345      IF (ierr .NE. NF_NOERR) THEN
346         PRINT*, "new_grid: Lecture echouee pour <rlonv>"
347         CALL abort
348      ENDIF
349c
350      ierr = NF_INQ_VARID (nid, "rlatu", nvarid)
351      IF (ierr .NE. NF_NOERR) THEN
352         PRINT*, "new_grid: Le champ <rlatu> est absent de start.nc"
353         CALL abort
354      ENDIF
355#ifdef NC_DOUBLE
356      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlatuold)
357#else
358      ierr = NF_GET_VAR_REAL(nid, nvarid, rlatuold)
359#endif
360      IF (ierr .NE. NF_NOERR) THEN
361         PRINT*, "new_grid: Lecture echouee pour <rlatu>"
362         CALL abort
363      ENDIF
364c
365      ierr = NF_INQ_VARID (nid, "rlonu", nvarid)
366      IF (ierr .NE. NF_NOERR) THEN
367         PRINT*, "new_grid: Le champ <rlonu> est absent de start.nc"
368         CALL abort
369      ENDIF
370#ifdef NC_DOUBLE
371      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlonuold)
372#else
373      ierr = NF_GET_VAR_REAL(nid, nvarid, rlonuold)
374#endif
375      IF (ierr .NE. NF_NOERR) THEN
376         PRINT*, "new_grid: Lecture echouee pour <rlonu>"
377         CALL abort
378      ENDIF
379c
380      ierr = NF_INQ_VARID (nid, "rlatv", nvarid)
381      IF (ierr .NE. NF_NOERR) THEN
382         PRINT*, "new_grid: Le champ <rlatv> est absent de start.nc"
383         CALL abort
384      ENDIF
385#ifdef NC_DOUBLE
386      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlatvold)
387#else
388      ierr = NF_GET_VAR_REAL(nid, nvarid, rlatvold)
389#endif
390      IF (ierr .NE. NF_NOERR) THEN
391         PRINT*, "new_grid: Lecture echouee pour <rlatv>"
392         CALL abort
393      ENDIF
394c
395
396      print*,"Lecture lon/lat OK"
397
398c-----------------------------------------------------------------------
399c Lecture des niveaux verticaux
400c-----------------------------------------------------------------------
401c
402      ierr = NF_INQ_VARID (nid, "nivsigs", nvarid)
403      IF (ierr .NE. NF_NOERR) THEN
404         PRINT*, "dynetat0: Le champ <nivsigs> est absent"
405         CALL abort
406      ENDIF
407#ifdef NC_DOUBLE
408      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, nivsigsold)
409#else
410      ierr = NF_GET_VAR_REAL(nid, nvarid, nivsigsold)
411#endif
412      IF (ierr .NE. NF_NOERR) THEN
413         PRINT*, "dynetat0: Lecture echouee pour <nivsigs>"
414         CALL abort
415      ENDIF
416
417      ierr = NF_INQ_VARID (nid, "nivsig", nvarid)
418      IF (ierr .NE. NF_NOERR) THEN
419         PRINT*, "dynetat0: Le champ <nivsig> est absent"
420         CALL abort
421      ENDIF
422#ifdef NC_DOUBLE
423      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, nivsigold)
424#else
425      ierr = NF_GET_VAR_REAL(nid, nvarid, nivsigold)
426#endif
427      IF (ierr .NE. NF_NOERR) THEN
428         PRINT*, "dynetat0: Lecture echouee pour <nivsig>"
429         CALL abort
430      ENDIF
431
432      ierr = NF_INQ_VARID (nid, "ap", nvarid)
433      IF (ierr .NE. NF_NOERR) THEN
434         PRINT*, "new_grid: Le champ <ap> est absent de start.nc"
435         CALL abort
436      ELSE
437#ifdef NC_DOUBLE
438         ierr = NF_GET_VAR_DOUBLE(nid, nvarid, apold)
439#else
440         ierr = NF_GET_VAR_REAL(nid, nvarid, apold)
441#endif
442         IF (ierr .NE. NF_NOERR) THEN
443            PRINT*, "new_grid: Lecture echouee pour <ap>"
444         ENDIF
445      ENDIF
446c
447      ierr = NF_INQ_VARID (nid, "bp", nvarid)
448      IF (ierr .NE. NF_NOERR) THEN
449         PRINT*, "new_grid: Le champ <bp> est absent de start.nc"
450         CALL abort
451      ENDIF
452#ifdef NC_DOUBLE
453      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, bpold)
454#else
455      ierr = NF_GET_VAR_REAL(nid, nvarid, bpold)
456#endif
457      IF (ierr .NE. NF_NOERR) THEN
458         PRINT*, "new_grid: Lecture echouee pour <bp>"
459         CALL abort
460      END IF
461
462      ierr = NF_INQ_VARID (nid, "presnivs", nvarid)
463      IF (ierr .NE. NF_NOERR) THEN
464         PRINT*, "dynetat0: Le champ <presnivs> est absent"
465         CALL abort
466      ENDIF
467#ifdef NC_DOUBLE
468      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, presnivsold)
469#else
470      ierr = NF_GET_VAR_REAL(nid, nvarid, presnivsold)
471#endif
472      IF (ierr .NE. NF_NOERR) THEN
473         PRINT*, "dynetat0: Lecture echouee pour <presnivs>"
474         CALL abort
475      ENDIF
476
477c-----------------------------------------------------------------------
478c Lecture geopotentiel au sol
479c-----------------------------------------------------------------------
480c
481      ierr = NF_INQ_VARID (nid, "phisinit", nvarid)
482      IF (ierr .NE. NF_NOERR) THEN
483         PRINT*, "new_grid: Le champ <phisinit> est absent de start.nc"
484         CALL abort
485      ENDIF
486#ifdef NC_DOUBLE
487      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, phisold)
488#else
489      ierr = NF_GET_VAR_REAL(nid, nvarid, phisold)
490#endif
491      IF (ierr .NE. NF_NOERR) THEN
492         PRINT*, "new_grid: Lecture echouee pour <phisinit>"
493         CALL abort
494      ENDIF
495
496      print*,"Lecture niveaux et geopot OK"
497
498c-----------------------------------------------------------------------
499c Lecture des champs 2D ()
500c-----------------------------------------------------------------------
501
502      start=(/1,1,1,0/)
503      counter=(/imold+1,jmold+1,1,0/)
504
505      ierr = NF_INQ_VARID (nid, "ps", nvarid)
506      IF (ierr .NE. NF_NOERR) THEN
507         PRINT*, "lect_start_archive: Le champ <ps> est absent"
508         CALL abort
509      ENDIF
510#ifdef NC_DOUBLE
511      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,psold)
512#else
513      ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,psold)
514#endif
515      IF (ierr .NE. NF_NOERR) THEN
516         PRINT*, "lect_start_archive: Lecture echouee pour <ps>"
517         CALL abort
518      ENDIF
519c
520      ierr = NF_INQ_VARID (nid, "tsurf", nvarid)
521      IF (ierr .NE. NF_NOERR) THEN
522         PRINT*, "lect_start_archive: Le champ <tsurf> est absent"
523         CALL abort
524      ENDIF
525#ifdef NC_DOUBLE
526      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,tsurfold)
527#else
528      ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,tsurfold)
529#endif
530      IF (ierr .NE. NF_NOERR) THEN
531         PRINT*, "lect_start_archive: Lecture echouee pour <tsurf>"
532         CALL abort
533      ENDIF
534c
535      do isoil=1,nsoilmx
536         write(str2,'(i2.2)') isoil
537c
538         ierr = NF_INQ_VARID (nid, "Tsoil"//str2, nvarid)
539         IF (ierr .NE. NF_NOERR) THEN
540            PRINT*, "lect_start_archive:
541     .              Le champ <","Tsoil"//str2,"> est absent"
542            CALL abort
543         ENDIF
544#ifdef NC_DOUBLE
545         ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,
546     .          tsoilold(1,1,isoil))
547#else
548         ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,
549     .          tsoilold(1,1,isoil))
550#endif
551         IF (ierr .NE. NF_NOERR) THEN
552            PRINT*, "lect_start_archive:
553     .            Lecture echouee pour <","Tsoil"//str2,">"
554            CALL abort
555         ENDIF
556      end do
557c
558      ierr = NF_INQ_VARID (nid, "albe", nvarid)
559      IF (ierr .NE. NF_NOERR) THEN
560         PRINT*, "lect_start_archive: Le champ <albe> est absent"
561         CALL abort
562      ENDIF
563#ifdef NC_DOUBLE
564      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,albeold)
565#else
566      ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,albeold)
567#endif
568      IF (ierr .NE. NF_NOERR) THEN
569         PRINT*, "lect_start_archive: Lecture echouee pour <albe>"
570         CALL abort
571      ENDIF
572c
573      ierr = NF_INQ_VARID (nid, "radsol", nvarid)
574      IF (ierr .NE. NF_NOERR) THEN
575         PRINT*, "lect_start_archive: Le champ <radsol> est absent"
576         CALL abort
577      ENDIF
578#ifdef NC_DOUBLE
579      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,radsolold)
580#else
581      ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,radsolold)
582#endif
583      IF (ierr .NE. NF_NOERR) THEN
584         PRINT*, "lect_start_archive: Lecture echouee pour <radsol>"
585         CALL abort
586      ENDIF
587c
588      ierr = NF_INQ_VARID (nid, "sollw", nvarid)
589      IF (ierr .NE. NF_NOERR) THEN
590         PRINT*, "lect_start_archive: Le champ <sollw> est absent"
591         CALL abort
592      ENDIF
593#ifdef NC_DOUBLE
594      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,sollwold)
595#else
596      ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,sollwold)
597#endif
598      IF (ierr .NE. NF_NOERR) THEN
599         PRINT*, "lect_start_archive: Lecture echouee pour <sollw>"
600         CALL abort
601      ENDIF
602c
603      ierr = NF_INQ_VARID (nid, "solsw", nvarid)
604      IF (ierr .NE. NF_NOERR) THEN
605         PRINT*, "lect_start_archive: Le champ <solsw> est absent"
606         CALL abort
607      ENDIF
608#ifdef NC_DOUBLE
609      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,solswold)
610#else
611      ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,solswold)
612#endif
613      IF (ierr .NE. NF_NOERR) THEN
614         PRINT*, "lect_start_archive: Lecture echouee pour <solsw>"
615         CALL abort
616      ENDIF
617c
618      ierr = NF_INQ_VARID (nid, "dlw", nvarid)
619      IF (ierr .NE. NF_NOERR) THEN
620         PRINT*, "lect_start_archive: Le champ <dlw> est absent"
621         CALL abort
622      ENDIF
623#ifdef NC_DOUBLE
624      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,dlwold)
625#else
626      ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,dlwold)
627#endif
628      IF (ierr .NE. NF_NOERR) THEN
629         PRINT*, "lect_start_archive: Lecture echouee pour <dlw>"
630         CALL abort
631      ENDIF
632c
633      ierr = NF_INQ_VARID (nid, "zmea", nvarid)
634      IF (ierr .NE. NF_NOERR) THEN
635         PRINT*, "lect_start_archive: Le champ <zmea> est absent"
636         PRINT*, "          Il est donc initialise a zero"
637         zmeaold=0.
638      ENDIF
639#ifdef NC_DOUBLE
640      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,zmeaold)
641#else
642      ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,zmeaold)
643#endif
644      IF (ierr .NE. NF_NOERR) THEN
645         PRINT*, "lect_start_archive: Lecture echouee pour <zmea>"
646         CALL abort
647      ENDIF
648c
649      ierr = NF_INQ_VARID (nid, "zstd", nvarid)
650      IF (ierr .NE. NF_NOERR) THEN
651         PRINT*, "lect_start_archive: Le champ <zstd> est absent"
652         PRINT*, "          Il est donc initialise a zero"
653         zstdold=0.
654      ENDIF
655#ifdef NC_DOUBLE
656      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,zstdold)
657#else
658      ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,zstdold)
659#endif
660      IF (ierr .NE. NF_NOERR) THEN
661         PRINT*, "lect_start_archive: Lecture echouee pour <zstd>"
662         CALL abort
663      ENDIF
664c
665      ierr = NF_INQ_VARID (nid, "zsig", nvarid)
666      IF (ierr .NE. NF_NOERR) THEN
667         PRINT*, "lect_start_archive: Le champ <zsig> est absent"
668         PRINT*, "          Il est donc initialise a zero"
669         zsigold=0.
670      ENDIF
671#ifdef NC_DOUBLE
672      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,zsigold)
673#else
674      ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,zsigold)
675#endif
676      IF (ierr .NE. NF_NOERR) THEN
677         PRINT*, "lect_start_archive: Lecture echouee pour <zsig>"
678         CALL abort
679      ENDIF
680c
681      ierr = NF_INQ_VARID (nid, "zgam", nvarid)
682      IF (ierr .NE. NF_NOERR) THEN
683         PRINT*, "lect_start_archive: Le champ <zgam> est absent"
684         PRINT*, "          Il est donc initialise a zero"
685         zgamold=0.
686      ENDIF
687#ifdef NC_DOUBLE
688      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,zgamold)
689#else
690      ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,zgamold)
691#endif
692      IF (ierr .NE. NF_NOERR) THEN
693         PRINT*, "lect_start_archive: Lecture echouee pour <zgam>"
694         CALL abort
695      ENDIF
696c
697      ierr = NF_INQ_VARID (nid, "zthe", nvarid)
698      IF (ierr .NE. NF_NOERR) THEN
699         PRINT*, "lect_start_archive: Le champ <zthe> est absent"
700         PRINT*, "          Il est donc initialise a zero"
701         ztheold=0.
702      ENDIF
703#ifdef NC_DOUBLE
704      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,ztheold)
705#else
706      ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,ztheold)
707#endif
708      IF (ierr .NE. NF_NOERR) THEN
709         PRINT*, "lect_start_archive: Lecture echouee pour <zthe>"
710         CALL abort
711      ENDIF
712c
713      ierr = NF_INQ_VARID (nid, "zpic", nvarid)
714      IF (ierr .NE. NF_NOERR) THEN
715         PRINT*, "lect_start_archive: Le champ <zpic> est absent"
716         PRINT*, "          Il est donc initialise a zero"
717         zpicold=0.
718      ENDIF
719#ifdef NC_DOUBLE
720      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,zpicold)
721#else
722      ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,zpicold)
723#endif
724      IF (ierr .NE. NF_NOERR) THEN
725         PRINT*, "lect_start_archive: Lecture echouee pour <zpic>"
726         CALL abort
727      ENDIF
728c
729      ierr = NF_INQ_VARID (nid, "zval", nvarid)
730      IF (ierr .NE. NF_NOERR) THEN
731         PRINT*, "lect_start_archive: Le champ <zval> est absent"
732         PRINT*, "          Il est donc initialise a zero"
733         zvalold=0.
734      ENDIF
735#ifdef NC_DOUBLE
736      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,zvalold)
737#else
738      ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,zvalold)
739#endif
740      IF (ierr .NE. NF_NOERR) THEN
741         PRINT*, "lect_start_archive: Lecture echouee pour <zval>"
742         CALL abort
743      ENDIF
744c
745
746      print*,"Lecture champs 2D OK"
747
748c-----------------------------------------------------------------------
749c       Lecture des champs 3D ()
750c-----------------------------------------------------------------------
751
752      start=(/1,1,1,1/)
753      counter=(/imold+1,jmold+1,lmold,1/)
754c
755      ierr = NF_INQ_VARID (nid,"temp", nvarid)
756      IF (ierr .NE. NF_NOERR) THEN
757         PRINT*, "lect_start_archive: Le champ <temp> est absent"
758         CALL abort
759      ENDIF
760#ifdef NC_DOUBLE
761      ierr = NF_GET_VARA_DOUBLE(nid, nvarid, start, counter, told)
762#else
763      ierr = NF_GET_VARA_REAL(nid, nvarid, start, counter, told)
764#endif
765      IF (ierr .NE. NF_NOERR) THEN
766         PRINT*, "lect_start_archive: Lecture echouee pour <temp>"
767         CALL abort
768      ENDIF
769c
770      ierr = NF_INQ_VARID (nid,"u", nvarid)
771      IF (ierr .NE. NF_NOERR) THEN
772         PRINT*, "lect_start_archive: Le champ <u> est absent"
773         CALL abort
774      ENDIF
775#ifdef NC_DOUBLE
776      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,uold)
777#else
778      ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,uold)
779#endif
780      IF (ierr .NE. NF_NOERR) THEN
781         PRINT*, "lect_start_archive: Lecture echouee pour <u>"
782         CALL abort
783      ENDIF
784c
785      ierr = NF_INQ_VARID (nid,"v", nvarid)
786      IF (ierr .NE. NF_NOERR) THEN
787         PRINT*, "lect_start_archive: Le champ <v> est absent"
788         CALL abort
789      ENDIF
790#ifdef NC_DOUBLE
791      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,vold)
792#else
793      ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,vold)
794#endif
795      IF (ierr .NE. NF_NOERR) THEN
796         PRINT*, "lect_start_archive: Lecture echouee pour <v>"
797         CALL abort
798      ENDIF
799c
800
801c TNAME: IL EST LU A PARTIR DE traceur.def (mettre l'ancien si
802c                changement du nombre de traceurs)
803
804      IF(iflag_trac.eq.1) THEN
805      DO iq=1,nqtot
806        ierr =  NF_INQ_VARID (nid, tname(iq), nvarid)
807        IF (ierr .NE. NF_NOERR) THEN
808            PRINT*, "dynetat0: Le champ <"//tname(iq)//"> est absent"
809            PRINT*, "          Il est donc initialise a zero"
810            qold(:,:,:,iq)=0.
811        ELSE
812#ifdef NC_DOUBLE
813      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,counter,qold(1,1,1,iq))
814#else
815      ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,qold(1,1,1,iq))
816#endif
817          IF (ierr .NE. NF_NOERR) THEN
818             PRINT*, "dynetat0: Lecture echouee pour "//tname(iq)
819             CALL abort
820          ENDIF
821        ENDIF
822      ENDDO
823      ENDIF
824
825
826      print*,"Lecture champs 3D OK"
827
828c=======================================================================
829c   INTERPOLATION DANS LA NOUVELLE GRILLE et initialisation des variables
830c=======================================================================
831c  Interpolation horizontale puis passage dans la grille physique pour
832c  les variables physique
833c  Interpolation verticale puis horizontale pour chaque variable 3D
834c=======================================================================
835
836c-----------------------------------------------------------------------
837c       Variable 2d :
838c-----------------------------------------------------------------------
839
840c Relief
841! topoflag = F: we keep the existing topography
842!          = T: we read it from the Relief.nc file
843! topoflag need to be in the specific run.def for newstart
844       topoflag = . FALSE .
845       CALL getin('topoflag',topoflag)
846
847       IF ( topoflag ) then
848         print*,"Topography (phis) read in file Relief.nc"
849         print*,"For Venus, map directly inverted in Relief.nc"
850         CALL startget_phys2d('surfgeo',iip1,jjp1,rlonv,rlatu,phis,0.0,
851     .                jjm ,rlonu,rlatv,.true.)
852c needed ?
853          phis(iip1,:) = phis(1,:)
854
855         CALL startget_phys1d('zmea',iip1,jjp1,rlonv,rlatu,ngridmx,zmea,
856     .            0.0,jjm,rlonu,rlatv,.true.)
857         CALL startget_phys1d('zstd',iip1,jjp1,rlonv,rlatu,ngridmx,zstd,
858     .            0.0,jjm,rlonu,rlatv,.true.)
859         CALL startget_phys1d('zsig',iip1,jjp1,rlonv,rlatu,ngridmx,zsig,
860     .            0.0,jjm,rlonu,rlatv,.true.)
861         CALL startget_phys1d('zgam',iip1,jjp1,rlonv,rlatu,ngridmx,zgam,
862     .            0.0,jjm,rlonu,rlatv,.true.)
863         CALL startget_phys1d('zthe',iip1,jjp1,rlonv,rlatu,ngridmx,zthe,
864     .            0.0,jjm,rlonu,rlatv,.true.)
865         CALL startget_phys1d('zpic',iip1,jjp1,rlonv,rlatu,ngridmx,zpic,
866     .            0.0,jjm,rlonu,rlatv,.true.)
867         CALL startget_phys1d('zval',iip1,jjp1,rlonv,rlatu,ngridmx,zval,
868     .            0.0,jjm,rlonu,rlatv,.true.)
869
870       ELSE
871          print*,'Using existing topography'
872          call interp_horiz (phisold,phis,imold,jmold,iim,jjm,1,
873     &                   rlonuold,rlatvold,rlonu,rlatv)
874
875          call interp_horiz (zmeaold,zmeaS,imold,jmold,iim,jjm,1,
876     &                   rlonuold,rlatvold,rlonu,rlatv)
877          call gr_dyn_fi (1,iip1,jjp1,ngridmx,zmeaS,zmea)
878          call interp_horiz (zstdold,zstdS,imold,jmold,iim,jjm,1,
879     &                   rlonuold,rlatvold,rlonu,rlatv)
880          call gr_dyn_fi (1,iip1,jjp1,ngridmx,zstdS,zstd)
881          call interp_horiz (zsigold,zsigS,imold,jmold,iim,jjm,1,
882     &                   rlonuold,rlatvold,rlonu,rlatv)
883          call gr_dyn_fi (1,iip1,jjp1,ngridmx,zsigS,zsig)
884          call interp_horiz (zgamold,zgamS,imold,jmold,iim,jjm,1,
885     &                   rlonuold,rlatvold,rlonu,rlatv)
886          call gr_dyn_fi (1,iip1,jjp1,ngridmx,zgamS,zgam)
887          call interp_horiz (ztheold,ztheS,imold,jmold,iim,jjm,1,
888     &                   rlonuold,rlatvold,rlonu,rlatv)
889          call gr_dyn_fi (1,iip1,jjp1,ngridmx,ztheS,zthe)
890          call interp_horiz (zpicold,zpicS,imold,jmold,iim,jjm,1,
891     &                   rlonuold,rlatvold,rlonu,rlatv)
892          call gr_dyn_fi (1,iip1,jjp1,ngridmx,zpicS,zpic)
893          call interp_horiz (zvalold,zvalS,imold,jmold,iim,jjm,1,
894     &                   rlonuold,rlatvold,rlonu,rlatv)
895          call gr_dyn_fi (1,iip1,jjp1,ngridmx,zvalS,zval)
896       ENDIF
897
898       print*,"New surface geopotential OK"
899
900c Lat et lon pour physique
901      do i=1,iip1
902        rlatS(i,:)=rlatu(:)*180./pi
903      enddo
904      call gr_dyn_fi (1,iip1,jjp1,ngridmx,rlatS,rlat)
905
906      do j=2,jjm
907        rlonS(:,j)=rlonv(:)*180./pi
908      enddo
909      rlonS(:,1)=0.
910      rlonS(:,jjp1)=0.
911      call gr_dyn_fi (1,iip1,jjp1,ngridmx,rlonS,rlon)
912
913c Temperature de surface
914      call interp_horiz (tsurfold,tsurfS,imold,jmold,iim,jjm,1,
915     &                   rlonuold,rlatvold,rlonu,rlatv)
916      call gr_dyn_fi (1,iip1,jjp1,ngridmx,tsurfS,tsurf)
917c     write(44,*) 'tsurf', tsurf
918
919c Temperature du sous-sol
920      call interp_horiz(tsoilold,tsoilS,
921     &                  imold,jmold,iim,jjm,nsoilmx,
922     &                   rlonuold,rlatvold,rlonu,rlatv)
923      call gr_dyn_fi (nsoilmx,iip1,jjp1,ngridmx,tsoilS,tsoil)
924c     write(45,*) 'tsoil',tsoil
925
926! CHANGING ALBEDO: may be done through run.def
927       albedoflag = . FALSE .
928       CALL getin('albedoflag',albedoflag)
929
930       IF ( albedoflag ) then
931         print*,"Albedo is reinitialized to the albedo value in run.def"
932         print*,"We may want to consider a map later on..."
933         albedo=0.1
934         CALL getin('albedo',albedo)
935         albe=albedo
936       ELSE
937         call interp_horiz (albeold,albeS,imold,jmold,iim,jjm,1,
938     &                   rlonuold,rlatvold,rlonu,rlatv)
939         call gr_dyn_fi (1,iip1,jjp1,ngridmx,albeS,albe)
940       ENDIF
941
942      call interp_horiz (radsolold,radsolS,imold,jmold,iim,jjm,1,
943     &                   rlonuold,rlatvold,rlonu,rlatv)
944      call gr_dyn_fi (1,iip1,jjp1,ngridmx,radsolS,radsol)
945
946      call interp_horiz (sollwold,sollwS,imold,jmold,iim,jjm,1,
947     &                   rlonuold,rlatvold,rlonu,rlatv)
948      call gr_dyn_fi (1,iip1,jjp1,ngridmx,sollwS,sollw)
949
950      call interp_horiz (solswold,solswS,imold,jmold,iim,jjm,1,
951     &                   rlonuold,rlatvold,rlonu,rlatv)
952      call gr_dyn_fi (1,iip1,jjp1,ngridmx,solswS,solsw)
953
954      call interp_horiz (dlwold,dlwS,imold,jmold,iim,jjm,1,
955     &                   rlonuold,rlatvold,rlonu,rlatv)
956      call gr_dyn_fi (1,iip1,jjp1,ngridmx,dlwS,dlw)
957
958      print*,"Nouvelles var physiques OK"
959
960c-----------------------------------------------------------------------
961c       Traitement special de la pression au sol :
962c-----------------------------------------------------------------------
963
964c  Extrapolation la pression dans la nouvelle grille
965      call interp_horiz(psold,ps,imold,jmold,iim,jjm,1,
966     &                   rlonuold,rlatvold,rlonu,rlatv)
967
968c       On assure la conservation de la masse de l'atmosphere
969c--------------------------------------------------------------
970
971      ptotal =  0.
972      DO j=1,jjp1
973         DO i=1,iim
974            ptotal=ptotal+ps(i,j)*aire(i,j)/g
975         END DO
976      END DO
977
978      write(*,*)
979      write(*,*)'Ancienne grille: masse de l atm :',ptotalold
980      write(*,*)'Nouvelle grille: masse de l atm :',ptotal
981      write (*,*) 'Ratio new atm./ old atm =', ptotal/ptotalold
982      write(*,*)
983
984
985      DO j=1,jjp1
986         DO i=1,iip1
987            ps(i,j)=ps(i,j) * ptotalold/ptotal
988         END DO
989      END DO
990
991c la pression de surface et les temperatures ne sont pas reequilibrees en fonction
992c de la nouvelle topographie...
993c Si l'ajustement inevitable du debut pose des problemes, voir le newstart martien.
994
995      print*,"Nouvelle ps OK"
996      print*,"If changes done on topography, beware !"
997      print*,"Some time may be needed for adjustments at the beginning"
998      print*,"so if unstable, relax the timestep and/or dissipation."
999
1000c-----------------------------------------------------------------------
1001c       Variable 3d :
1002c-----------------------------------------------------------------------
1003
1004      CALL pression(ip1jmp1, ap, bp, ps, p3d)
1005         if (disvert_type==1) then
1006           CALL exner_hyb(  ip1jmp1, ps, p3d,alpha,beta,pks, pk, pkf )
1007         else ! we assume that we are in the disvert_type==2 case
1008           CALL exner_milieu( ip1jmp1, ps, p3d, beta, pks, pk, pkf )
1009         endif
1010     
1011c temperatures atmospheriques
1012
1013c ATTENTION: peut servir, mais bon...
1014c modif: profil uniforme
1015c     do l=1,lmold
1016c      do j=1,jmold+1
1017c       do i=1,imold+1
1018c          told(i,j,l)=told(1,jmold/2,l)
1019c       enddo
1020c      enddo
1021c     enddo
1022
1023      write (*,*) 'told ', told (1,jmold+1,1)  ! INFO
1024      call interp_vert
1025     &    (told,var,lmold,llm,apold,bpold,ap,bp,
1026     &     psold,(imold+1)*(jmold+1))
1027      write (*,*) 'var ', var (1,jmold+1,1)  ! INFO
1028      call interp_horiz(var,t,imold,jmold,iim,jjm,llm,
1029     &                   rlonuold,rlatvold,rlonu,rlatv)
1030      write (*,*) 'T ', T(1,jjp1,1)  ! INFO
1031! pour info:
1032! Si extension verticale, la T est extrapolee constante au-dessus de lmold
1033
1034c passage grille physique pour restartphy.nc
1035      call gr_dyn_fi (llm,iip1,jjp1,ngridmx,T,t_fi)
1036
1037! ADAPTATION GCM POUR CP(T)
1038c passage en temperature potentielle
1039      call t2tpot(ip1jmp1*llm,T,teta,pk)
1040c on assure la periodicite
1041      teta(iip1,:,:) =  teta(1,:,:)
1042
1043c calcul des champ de vent; passage en vent covariant
1044      write (*,*) 'uold ', uold (1,2,1)  ! INFO
1045      call interp_vert
1046     & (uold,var,lmold,llm,apold,bpold,ap,bp,
1047     &  psold,(imold+1)*(jmold+1))
1048      write (*,*) 'var ', var (1,2,1)  ! INFO
1049      call interp_horiz(var,us,imold,jmold,iim,jjm,llm,
1050     &                   rlonuold,rlatvold,rlonu,rlatv)
1051      write (*,*) 'us ', us (1,2,1)   ! INFO
1052
1053      call interp_vert
1054     & (vold,var,lmold,llm,
1055     &  apold,bpold,ap,bp,psold,(imold+1)*(jmold+1))
1056      call interp_horiz(var,vs,imold,jmold,iim,jjm,llm,
1057     &                   rlonuold,rlatvold,rlonu,rlatv)
1058      call scal_wind(us,vs,unat,vnat)
1059      write (*,*) 'unat ', unat (1,2,1)    ! INFO
1060      do l=1,llm
1061        do j = 1, jjp1
1062          do i=1,iip1
1063            ucov( i,j,l ) = unat( i,j,l ) * cu(i,j)
1064! pour info:
1065! Si extension verticale, on impose u=0 au-dessus de lmold
1066            if (l.gt.lmold) ucov( i,j,l ) = 0
1067          end do
1068        end do
1069      end do
1070      write (*,*) 'ucov ', ucov (1,2,1)  ! INFO
1071c     write(48,*) 'ucov',ucov
1072      do l=1,llm
1073        do j = 1, jjm
1074          do i=1,iim
1075            vcov( i,j,l ) = vnat( i,j,l ) * cv(i,j)
1076! pour info:
1077! Si extension verticale, on impose v=0 au-dessus de lmold
1078            if (l.gt.lmold) vcov( i,j,l ) = 0
1079          end do
1080          vcov( iip1,j,l ) = vcov( 1,j,l )
1081        end do
1082      end do
1083c     write(49,*) 'ucov',vcov
1084
1085c masse: on la recalcule (ps a été ajustée pour conserver la masse totale)
1086      call massdair(p3d,masse)
1087     
1088c traceurs 3D
1089      do  iq = 1, nqtot
1090            call interp_vert(qold(1,1,1,iq),var,lmold,llm,
1091     &        apold,bpold,ap,bp,psold,(imold+1)*(jmold+1))
1092            call interp_horiz(var,q(1,1,1,iq),imold,jmold,iim,jjm,llm,
1093     &                  rlonuold,rlatvold,rlonu,rlatv)
1094      enddo
1095
1096      print*,"Nouvelles var dynamiques OK"
1097
1098c=======================================================================
1099c    Ecriture des restart.nc et restartphy.nc :
1100c=======================================================================
1101
1102      call writerestart('restart.nc',tab_cntrl_dyn,
1103     .                    phis,vcov,ucov,teta,masse,q,ps)
1104
1105      print*,"restart OK"
1106
1107      call writerestartphy('restartphy.nc',tab_cntrl_fi,ngridmx,llm,
1108     .           rlat,rlon,tsurf,tsoil,albe,solsw, sollw,dlw,
1109     .           radsol,
1110     .           zmea,zstd,zsig,zgam,zthe,zpic,zval,
1111     .           t_fi)
1112
1113      print*,"restartphy OK"
1114
1115      end
Note: See TracBrowser for help on using the repository browser.