source: trunk/LMDZ.COMMON/libf/dyn3d/newstart-VT.F @ 825

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