source: trunk/LMDZ.COMMON/libf/dynlonlat_phylonlat/phytitan/newstart.F @ 1428

Last change on this file since 1428 was 1403, checked in by emillour, 11 years ago

All models: Reorganizing the physics/dynamics interface.

  • makelmdz and makelmdz_fcm scripts adapted to handle the new directory settings
  • misc: (replaces what was the "bibio" directory)
  • Should only contain extremely generic (and non physics or dynamics-specific) routines
  • Therefore moved initdynav.F90, initfluxsto.F, inithist.F, writedynav.F90, write_field.F90, writehist.F to "dyn3d_common"
  • dynlonlat_phylonlat: (new interface directory)
  • This directory contains routines relevent to physics/dynamics grid interactions, e.g. routines gr_dyn_fi or gr_fi_dyn and calfis
  • Moreover the dynlonlat_phylonlat contains directories "phy*" corresponding to each physics package "phy*" to be used. These subdirectories should only contain specific interfaces (e.g. iniphysiq) or main programs (e.g. newstart)
  • phy*/dyn1d: this subdirectory contains the 1D model using physics from phy*

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