source: trunk/LMDZ.COMMON/libf/dynlonlat_phylonlat/phyvenus/newstart.F @ 1455

Last change on this file since 1455 was 1443, checked in by emillour, 10 years ago

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