source: trunk/mars/libf/dyn3d/lect_start_archive.F @ 77

Last change on this file since 77 was 38, checked in by emillour, 14 years ago

Ajout du modè Martien (mon LMDZ.MARS.BETA, du 28/01/2011) dans le rértoire mars, pour pouvoir suivre plus facilement les modifs.
EM
File size: 48.9 KB
RevLine 
[38]1      SUBROUTINE lect_start_archive(date,tsurf,tsoil,emis,q2,
2     &     t,ucov,vcov,ps,co2ice,h,phisold_newgrid,
3     &     q,qsurf,surfith,nid)
4c=======================================================================
5c
6c
7c   Auteur:    05/1997 , 12/2003 : coord hybride  FF
8c   ------
9c
10c
11c   Objet:     Lecture des variables d'un fichier "start_archive"
12c              Plus besoin de régler ancienne valeurs grace
13c              a l'allocation dynamique de memoire (Yann Wanherdrick)
14c
15c
16c
17c=======================================================================
18
19      implicit none
20
21#include "dimensions.h"
22#include "dimphys.h"
23#include "surfdat.h"
24#include "comsoil.h"
25#include "dimradmars.h"
26#include "yomaer.h"
27#include "planete.h"
28#include "paramet.h"
29#include "comconst.h"
30#include "comvert.h"
31#include "comgeom2.h"
32#include "control.h"
33#include "logic.h"
34#include "description.h"
35#include "ener.h"
36#include "temps.h"
37#include "lmdstd.h"
38#include "netcdf.inc"
39!#include "tracer.h"
40#include"advtrac.h"
41c=======================================================================
42c   Declarations
43c=======================================================================
44
45c Old variables dimensions (from file)
46c------------------------------------
47      INTEGER   imold,jmold,lmold,nsoilold,nqold
48
49c et autres:
50c----------
51      INTEGER lnblnk
52      EXTERNAL lnblnk
53
54c Variables pour les lectures des fichiers "ini"
55c--------------------------------------------------
56!      INTEGER sizei,
57      integer timelen,dimid
58!      INTEGER length
59!      parameter (length = 100)
60      INTEGER tab0
61      INTEGER isoil,iq,iqmax
62      CHARACTER*2   str2
63
64!      REAL dimfirst(4) ! tableau contenant les 1ers elements des dimensions
65
66!      REAL dimlast(4) ! tableau contenant les derniers elements des dimensions
67
68!      REAL dimcycl(4) ! tableau contenant les periodes des dimensions
69!      CHARACTER*120 dimsource
70!      CHARACTER*16 dimname
71!      CHARACTER*80 dimtitle
72!      CHARACTER*40 dimunits
73!      INTEGER   dimtype
74
75!      INTEGER dimord(4)  ! tableau contenant l''ordre
76!      data dimord /1,2,3,4/ ! de sortie des dimensions
77
78!      INTEGER vardim(4)
79      REAL date
80      INTEGER   memo
81!      character (len=50) :: tmpname
82
83c Variable histoire
84c------------------
85      REAL vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm) ! vents covariants
86      REAL h(iip1,jjp1,llm),ps(iip1,jjp1)
87      REAL q(iip1,jjp1,llm,nqmx),qtot(iip1,jjp1,llm)
88
89c autre variables dynamique nouvelle grille
90c------------------------------------------
91
92c!-*-
93!      integer klatdat,klongdat
94!      PARAMETER (klatdat=180,klongdat=360)
95
96c Physique sur grille scalaire
97c----------------------------
98
99c variable physique
100c------------------
101      REAL tsurf(ngridmx) ! surface temperature
102      REAL tsoil(ngridmx,nsoilmx) ! soil temperature
103      REAL co2ice(ngridmx) ! CO2 ice layer
104      REAL emis(ngridmx)
105      REAL q2(ngridmx,nlayermx+1),qsurf(ngridmx,nqmx)
106c     REAL phisfi(ngridmx)
107
108      INTEGER i,j,l
109      INTEGER nid,nvarid
110c     REAL year_day,periheli,aphelie,peri_day
111c     REAL obliquit,z0,emin_turb,lmixmin
112c     REAL emissiv,emisice(2),albedice(2),tauvis
113c     REAL iceradius(2) , dtemisice(2)
114
115!      EXTERNAL RAN1
116!      REAL RAN1
117!      EXTERNAL geopot,inigeom
118      integer ierr
119!      integer ismin
120!      external ismin
121!      CHARACTER*80 datapath
122      integer, dimension(4) :: start,count
123
124c Variable nouvelle grille naturelle au point scalaire
125c------------------------------------------------------
126      real us(iip1,jjp1,llm),vs(iip1,jjp1,llm)
127      REAL phisold_newgrid(iip1,jjp1)
128      REAL t(iip1,jjp1,llm)
129      real tsurfS(iip1,jjp1),tsoilS(iip1,jjp1,nsoilmx)
130      real inertiedatS(iip1,jjp1,nsoilmx)
131      real co2iceS(iip1,jjp1),emisS(iip1,jjp1)
132      REAL q2S(iip1,jjp1,llm+1),qsurfS(iip1,jjp1,nqmx)
133
134      real ptotal, co2icetotal
135
136c Var intermediaires : vent naturel, mais pas coord scalaire
137c-----------------------------------------------------------
138      real vnat(iip1,jjm,llm),unat(iip1,jjp1,llm)
139
140
141c Variable de l'ancienne grille
142c---------------------------------------------------------
143
144      real, dimension(:), allocatable :: timelist
145      real, dimension(:), allocatable :: rlonuold, rlatvold
146      real, dimension(:), allocatable :: rlonvold, rlatuold
147      real, dimension(:), allocatable :: apsold,bpsold
148      real, dimension(:), allocatable :: mlayerold
149      real, dimension(:,:,:), allocatable :: uold,vold,told,q2old
150      real, dimension(:,:,:), allocatable :: tsoilold,qsurfold
151      real, dimension(:,:,:),allocatable :: tsoiloldnew
152! tsoiloldnew: old soil values, but along new subterranean grid
153      real, dimension(:,:,:), allocatable :: inertiedatold
154! inertiedatoldnew: old inertia values, but along new subterranean grid
155      real, dimension(:,:,:), allocatable :: inertiedatoldnew
156      real, dimension(:,:), allocatable :: psold,phisold
157      real, dimension(:,:), allocatable :: co2iceold,tsurfold
158      real, dimension(:,:), allocatable :: emisold
159      real, dimension(:,:,:,:), allocatable :: qold
160
161      real tab_cntrl(100)
162
163      real ptotalold, co2icetotalold
164
165      logical :: olddepthdef=.false. ! flag
166! olddepthdef=.true. if soil depths are in 'old' (unspecified) format
167      logical :: depthinterpol=.false. ! flag
168! depthinterpol=.true. if interpolation will be requiered
169      logical :: therminertia_3D=.true. ! flag
170! therminertia_3D=.true. if thermal inertia is 3D and read from datafile
171c Variable intermediaires iutilise pour l'extrapolation verticale
172c----------------------------------------------------------------
173      real, dimension(:,:,:), allocatable :: var,varp1
174      real, dimension(:), allocatable :: oldgrid, oldval
175      real, dimension(:), allocatable :: newval
176!      real, dimension(:), allocatable :: oldmlayer
177
178      real surfith(iip1,jjp1) ! surface thermal inertia
179!      real surfithfi(ngridmx)
180      ! surface thermal inertia at old horizontal grid resolution
181      real, dimension(:,:), allocatable :: surfithold
182
183! flag which identifies if archive file is using old tracer names (qsurf01,...)
184      logical :: oldtracernames=.false.
185      integer :: counter
186      character(len=30) :: txt ! to store some text
187
188c=======================================================================
189
190! 0. Preliminary stuff
191
192! check if tracers follow old naming convention (q01, q02, q03, ...)
193      counter=0
194      do iq=1,nqmx
195        txt= " "
196        write(txt,'(a1,i2.2)')'q',iq
197        ierr=NF_INQ_VARID(nid,txt,nvarid)
198        if (ierr.ne.NF_NOERR) then
199          ! did not find old tracer name
200          exit ! might as well stop here
201        else
202          ! found old tracer name
203          counter=counter+1
204        endif
205      enddo
206      if (counter.eq.nqmx) then
207        write(*,*) "lect_start_archive: tracers seem to follow old ",
208     &             "naming convention (q01, q02,...)"
209        oldtracernames=.true.
210      endif
211
212
213!-----------------------------------------------------------------------
214! 1. Read data dimensions (i.e. size and length)
215!-----------------------------------------------------------------------
216
217! 1.2 Read the various dimension lengths of data in file
218
219      ierr= NF_INQ_DIMID(nid,"Time",dimid)
220      if (ierr.ne.NF_NOERR) then
221         ierr= NF_INQ_DIMID(nid,"temps",dimid)
222      endif
223      ierr= NF_INQ_DIMLEN(nid,dimid,timelen)
224
225      ierr= NF_INQ_DIMID(nid,"latitude",dimid)
226      if (ierr.ne.NF_NOERR) then
227         ierr= NF_INQ_DIMID(nid,"rlatu",dimid)
228      endif
229      ierr= NF_INQ_DIMLEN(nid,dimid,jmold)
230      jmold=jmold-1
231
232      ierr= NF_INQ_DIMID(nid,"longitude",dimid)
233      if (ierr.ne.NF_NOERR) then
234         ierr= NF_INQ_DIMID(nid,"rlonv",dimid)
235      endif
236      ierr= NF_INQ_DIMLEN(nid,dimid,imold)
237      imold=imold-1
238
239      ierr= NF_INQ_DIMID(nid,"altitude",dimid)
240      if (ierr.ne.NF_NOERR) then
241         ierr= NF_INQ_DIMID(nid,"sig_s",dimid)
242      endif
243      ierr= NF_INQ_DIMLEN(nid,dimid,lmold)
244
245      nqold=0
246      do
247         write(str2,'(i2.2)') nqold+1
248         ierr= NF_INQ_VARID(nid,'q'//str2,dimid)
249!        write(*,*) 'q'//str2
250         if (ierr.eq.NF_NOERR) then
251            nqold=nqold+1
252         else
253            exit
254         endif
255      enddo
256
257! 1.2.1 find out the # of subsurface_layers
258      nsoilold=0 !dummy initialisation
259      ierr=NF_INQ_DIMID(nid,"subsurface_layers",dimid)
260      if (ierr.eq.NF_NOERR) then
261        ierr=NF_INQ_DIMLEN(nid,dimid,nsoilold)
262        if (ierr.ne.NF_NOERR) then
263         write(*,*)'lec_start_archive: ',
264     &              'Failed reading subsurface_layers length'
265        endif
266      else
267        write(*,*)"lec_start_archive: did not find subsurface_layers"
268      endif
269
270      if (nsoilold.eq.0) then ! 'old' archive format;
271      ! must use Tg//str2 fields to compute nsoilold
272      write(*,*)"lec_start_archive: building nsoilold from Tg fields"
273        do
274         write(str2,'(i2.2)') nsoilold+1
275         ierr=NF_INQ_VARID(nid,'Tg'//str2,dimid)
276         if (ierr.eq.NF_NOERR) then
277          nsoilold=nsoilold+1
278         else
279          exit
280         endif
281        enddo
282      endif
283
284
285      if (nsoilold.ne.nsoilmx) then ! interpolation will be required
286        depthinterpol=.true.
287      endif
288
289! 1.3 Report dimensions
290     
291      write(*,*) "Start_archive dimensions:"
292      write(*,*) "longitude: ",imold
293      write(*,*) "latitude: ",jmold
294      write(*,*) "altitude: ",lmold
295      write(*,*) "tracers: ",nqold
296      write(*,*) "subsurface_layers: ",nsoilold
297      if (depthinterpol) then
298      write(*,*) " => Warning, nsoilmx= ",nsoilmx
299      write(*,*) '    which implies that you want subterranean interpola
300     &tion.'
301      write(*,*) '  Otherwise, set nsoilmx -in dimphys.h- to: ',nsoilold
302      endif
303      write(*,*) "time lenght: ",timelen
304      write(*,*)
305
306!-----------------------------------------------------------------------
307! 2. Allocate arrays to store datasets
308!-----------------------------------------------------------------------
309
310      allocate(timelist(timelen))
311      allocate(rlonuold(imold+1), rlatvold(jmold))
312      allocate(rlonvold(imold+1), rlatuold(jmold+1))
313      allocate (apsold(lmold),bpsold(lmold))
314      allocate(uold(imold+1,jmold+1,lmold))
315      allocate(vold(imold+1,jmold+1,lmold))
316      allocate(told(imold+1,jmold+1,lmold))
317      allocate(psold(imold+1,jmold+1))
318      allocate(phisold(imold+1,jmold+1))
319      allocate(qold(imold+1,jmold+1,lmold,nqmx))
320      allocate(co2iceold(imold+1,jmold+1))
321      allocate(tsurfold(imold+1,jmold+1))
322      allocate(emisold(imold+1,jmold+1))
323      allocate(q2old(imold+1,jmold+1,lmold+1))
324!      allocate(tsoilold(imold+1,jmold+1,nsoilmx))
325      allocate(tsoilold(imold+1,jmold+1,nsoilold))
326      allocate(tsoiloldnew(imold+1,jmold+1,nsoilmx))
327      allocate(inertiedatold(imold+1,jmold+1,nsoilold)) ! soil thermal inertia
328      allocate(inertiedatoldnew(imold+1,jmold+1,nsoilmx))
329      ! surface thermal inertia at old horizontal grid resolution
330      allocate(surfithold(imold+1,jmold+1))
331      allocate(mlayerold(nsoilold))
332      allocate(qsurfold(imold+1,jmold+1,nqmx))
333
334      allocate(var (imold+1,jmold+1,llm))
335      allocate(varp1 (imold+1,jmold+1,llm+1))
336
337      write(*,*) 'q2',ngridmx,nlayermx+1
338      write(*,*) 'q2S',iip1,jjp1,llm+1
339      write(*,*) 'q2old',imold+1,jmold+1,lmold+1
340
341!-----------------------------------------------------------------------
342! 3. Read time-independent data
343!-----------------------------------------------------------------------
344
345C-----------------------------------------------------------------------
346c 3.1. Lecture du tableau des parametres du run
347c     (pour  la lecture ulterieure de "ptotalold" et "co2icetotalold")
348c-----------------------------------------------------------------------
349c
350      ierr = NF_INQ_VARID (nid, "controle", nvarid)
351      IF (ierr .NE. NF_NOERR) THEN
352         PRINT*, "Lect_start_archive: champ <controle> est absent"
353         CALL abort
354      ENDIF
355#ifdef NC_DOUBLE
356      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tab_cntrl)
357#else
358      ierr = NF_GET_VAR_REAL(nid, nvarid, tab_cntrl)
359#endif
360      IF (ierr .NE. NF_NOERR) THEN
361         PRINT*, "lect_start_archive: Lecture echoue pour <controle>"
362         CALL abort
363      ENDIF
364c
365      tab0 = 50
366
367c-----------------------------------------------------------------------
368c 3.2 Lecture des longitudes et latitudes
369c-----------------------------------------------------------------------
370c
371      ierr = NF_INQ_VARID (nid, "rlonv", nvarid)
372      IF (ierr .NE. NF_NOERR) THEN
373         PRINT*, "lect_start_archive: Le champ <rlonv> est absent"
374         CALL abort
375      ENDIF
376#ifdef NC_DOUBLE
377      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlonvold)
378#else
379      ierr = NF_GET_VAR_REAL(nid, nvarid, rlonvold)
380#endif
381      IF (ierr .NE. NF_NOERR) THEN
382         PRINT*, "lect_start_archive: Lecture echouee pour <rlonv>"
383         CALL abort
384      ENDIF
385c
386      ierr = NF_INQ_VARID (nid, "rlatu", nvarid)
387      IF (ierr .NE. NF_NOERR) THEN
388         PRINT*, "lect_start_archive: Le champ <rlatu> est absent"
389         CALL abort
390      ENDIF
391#ifdef NC_DOUBLE
392      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlatuold)
393#else
394      ierr = NF_GET_VAR_REAL(nid, nvarid, rlatuold)
395#endif
396      IF (ierr .NE. NF_NOERR) THEN
397         PRINT*, "lect_start_archive: Lecture echouee pour <rlatu>"
398         CALL abort
399      ENDIF
400c
401      ierr = NF_INQ_VARID (nid, "rlonu", nvarid)
402      IF (ierr .NE. NF_NOERR) THEN
403         PRINT*, "lect_start_archive: Le champ <rlonu> est absent"
404         CALL abort
405      ENDIF
406#ifdef NC_DOUBLE
407      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlonuold)
408#else
409      ierr = NF_GET_VAR_REAL(nid, nvarid, rlonuold)
410#endif
411      IF (ierr .NE. NF_NOERR) THEN
412         PRINT*, "lect_start_archive: Lecture echouee pour <rlonu>"
413         CALL abort
414      ENDIF
415c
416      ierr = NF_INQ_VARID (nid, "rlatv", nvarid)
417      IF (ierr .NE. NF_NOERR) THEN
418         PRINT*, "lect_start_archive: Le champ <rlatv> est absent"
419         CALL abort
420      ENDIF
421#ifdef NC_DOUBLE
422      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlatvold)
423#else
424      ierr = NF_GET_VAR_REAL(nid, nvarid, rlatvold)
425#endif
426      IF (ierr .NE. NF_NOERR) THEN
427         PRINT*, "lect_start_archive: Lecture echouee pour <rlatv>"
428         CALL abort
429      ENDIF
430c
431
432c-----------------------------------------------------------------------
433c 3.3. Lecture des niveaux verticaux
434c-----------------------------------------------------------------------
435c
436      ierr = NF_INQ_VARID (nid, "aps", nvarid)
437      IF (ierr .NE. NF_NOERR) THEN
438         PRINT*, "lect_start_archive: Le champ <aps> est absent"
439         apsold=0
440         PRINT*, "<aps> set to 0"
441      ELSE
442#ifdef NC_DOUBLE
443         ierr = NF_GET_VAR_DOUBLE(nid, nvarid, apsold)
444#else
445         ierr = NF_GET_VAR_REAL(nid, nvarid, apsold)
446#endif
447         IF (ierr .NE. NF_NOERR) THEN
448            PRINT*, "lect_start_archive: Lecture echouee pour <aps>"
449         ENDIF
450      ENDIF
451c
452      ierr = NF_INQ_VARID (nid, "bps", nvarid)
453      IF (ierr .NE. NF_NOERR) THEN
454         PRINT*, "lect_start_archive: Le champ <bps> est absent"
455         PRINT*, "It must be an old start_archive, lets look for sig_s"
456         ierr = NF_INQ_VARID (nid, "sig_s", nvarid)
457         IF (ierr .NE. NF_NOERR) THEN
458            PRINT*, "Nothing to do..."
459            CALL abort
460         ENDIF
461      ENDIF
462#ifdef NC_DOUBLE
463      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, bpsold)
464#else
465      ierr = NF_GET_VAR_REAL(nid, nvarid, bpsold)
466#endif
467      IF (ierr .NE. NF_NOERR) THEN
468         PRINT*, "lect_start_archive: Lecture echouee pour <bps>"
469         CALL abort
470      END IF
471
472c-----------------------------------------------------------------------
473c 3.4 Read Soil layers depths
474c-----------------------------------------------------------------------
475     
476      ierr=NF_INQ_VARID(nid,"soildepth",nvarid)
477      if (ierr.ne.NF_NOERR) then
478       write(*,*)'lect_start_archive: Could not find <soildepth>'
479       write(*,*)' => Assuming this is an archive in old format'
480       olddepthdef=.true.
481       depthinterpol=.true.
482       ! this is how soil depth was defined in ye old days
483        do isoil=1,nsoilold
484          mlayerold(isoil)=sqrt(887.75/3.14)*((2.**(isoil-0.5))-1.)
485        enddo
486      else
487#ifdef NC_DOUBLE
488        ierr = NF_GET_VAR_DOUBLE(nid,nvarid,mlayerold)
489#else
490        ierr = NF_GET_VAR_REAL(nid,nvarid,mlayerold)
491#endif
492       if (ierr .NE. NF_NOERR) then
493         PRINT*, "lect_start_archive: Failed reading <soildepth>"
494         CALL abort
495       endif
496
497      endif !of if(ierr.ne.NF_NOERR)
498
499      ! Read (or build) mlayer()
500      if (depthinterpol) then
501       ! Build (default) new soil depths (mlayer(:) is in comsoil.h),
502       ! as in soil_settings.F
503       write(*,*)' => Building default soil depths'
504       do isoil=0,nsoilmx-1
505         mlayer(isoil)=2.e-4*(2.**(isoil-0.5))
506       enddo
507       write(*,*)' => mlayer: ',mlayer
508       ! Also build (default) new soil interlayer depth layer(:)
509       do isoil=1,nsoilmx
510         layer(isoil)=sqrt(mlayer(0)*mlayer(1))*
511     &                      ((mlayer(1)/mlayer(0))**(isoil-1))
512       enddo
513       write(*,*)' =>  layer: ',layer
514      else ! read mlayer() from file
515#ifdef NC_DOUBLE
516        ierr = NF_GET_VAR_DOUBLE(nid,nvarid,mlayer)
517#else
518        ierr = NF_GET_VAR_REAL(nid,nvarid,mlayer)
519#endif
520       if (ierr .NE. NF_NOERR) then
521         PRINT*, "lect_start_archive: Failed reading <soildepth>"
522         CALL abort
523       endif
524      endif ! of if (depthinterpol)
525
526c-----------------------------------------------------------------------
527c 3.5 Read Soil thermal inertia
528c-----------------------------------------------------------------------
529
530      ierr=NF_INQ_VARID(nid,"inertiedat",nvarid)
531      if (ierr.ne.NF_NOERR) then
532       write(*,*)'lect_start_archive: Could not find <inertiedat>'
533       write(*,*)' => Assuming this is an archive in old format'
534       therminertia_3D=.false.
535       write(*,*)' => Thermal inertia will be read from reference file'
536       volcapa=1.e6
537       write(*,*)'    and soil volumetric heat capacity is set to ',
538     &           volcapa
539      else
540#ifdef NC_DOUBLE
541        ierr = NF_GET_VAR_DOUBLE(nid,nvarid,inertiedatold)
542#else
543        ierr = NF_GET_VAR_REAL(nid,nvarid,inertiedatold)
544#endif
545       if (ierr .NE. NF_NOERR) then
546         PRINT*, "lect_start_archive: Failed reading <inertiedat>"
547         CALL abort
548       endif
549      endif
550
551c-----------------------------------------------------------------------
552c 3.6 Lecture geopotentiel au sol
553c-----------------------------------------------------------------------
554c
555      ierr = NF_INQ_VARID (nid, "phisinit", nvarid)
556      IF (ierr .NE. NF_NOERR) THEN
557         PRINT*, "lect_start_archive: Le champ <phisinit> est absent"
558         CALL abort
559      ENDIF
560#ifdef NC_DOUBLE
561      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, phisold)
562#else
563      ierr = NF_GET_VAR_REAL(nid, nvarid, phisold)
564#endif
565      IF (ierr .NE. NF_NOERR) THEN
566         PRINT*, "lect_start_archive: Lecture echouee pour <phisinit>"
567         CALL abort
568      ENDIF
569
570C-----------------------------------------------------------------------
571c   lecture de "ptotalold" et "co2icetotalold"
572c-----------------------------------------------------------------------
573      ptotalold = tab_cntrl(tab0+49)
574      co2icetotalold = tab_cntrl(tab0+50)
575 
576c-----------------------------------------------------------------------
577c 4. Lecture du temps et choix
578c-----------------------------------------------------------------------
579 
580c  lecture du temps
581c
582      ierr = NF_INQ_DIMID (nid, "Time", nvarid)
583      IF (ierr .NE. NF_NOERR) THEN
584         ierr = NF_INQ_DIMID (nid, "temps", nvarid)
585         IF (ierr .NE. NF_NOERR) THEN
586            PRINT*, "lect_start_archive: Le champ <Time> est absent"
587            CALL abort
588         endif
589      ENDIF
590
591      ierr = NF_INQ_VARID (nid, "Time", nvarid)
592      IF (ierr .NE. NF_NOERR) THEN
593         ierr = NF_INQ_VARID (nid, "temps", nvarid)
594      endif
595#ifdef NC_DOUBLE
596      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, timelist)
597#else
598      ierr = NF_GET_VAR_REAL(nid, nvarid, timelist)
599#endif
600      IF (ierr .NE. NF_NOERR) THEN
601         PRINT*, "lect_start_archive: Lecture echouee pour <Time>"
602         CALL abort
603      ENDIF
604c
605      write(*,*)
606      write(*,*)
607      write(*,*) 'Differentes dates des etats initiaux stockes:'
608      write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
609      pi=2.*ASIN(1.)
610      do i=1,timelen
611c       call solarlong(timelist(i),sollong(i))
612c       sollong(i) = sollong(i)*180./pi
613        write(*,*) 'etat initial au jour martien' ,int(timelist(i))
614c       write(*,6) nint(timelist(i)),nint(mod(timelist(i),669)),
615c    .    sollong(i)
616      end do
617
618   6  FORMAT(i7,i7,f9.3)
619 
620      write(*,*)
621      write(*,*) 'Choix de la date'
622 123  read(*,*,iostat=ierr) date
623      if(ierr.ne.0) goto 123
624      memo = 0
625      do i=1,timelen
626        if (date.eq.int(timelist(i))) then
627            memo = i
628        endif
629      end do
630 
631      if (memo.eq.0) then
632        write(*,*)
633        write(*,*)
634        write(*,*) 'He alors... Y sait pas lire !?!'
635        write(*,*)
636        write(*,*) 'Differentes dates des etats initiaux stockes:'
637        write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
638        do i=1,timelen
639          write(*,*) 'etat initial au jour martien' ,nint(timelist(i))
640c         write(*,6) nint(timelist(i)),nint(mod(timelist(i),669))
641        end do
642        goto 123
643      endif
644
645!-----------------------------------------------------------------------
646! 5. Read (time-dependent) data from datafile
647!-----------------------------------------------------------------------
648
649
650c-----------------------------------------------------------------------
651c 5.1 Lecture des champs 2D (co2ice, emis,ps,tsurf,Tg[10], q2surf)
652c-----------------------------------------------------------------------
653 
654      start=(/1,1,memo,0/)
655      count=(/imold+1,jmold+1,1,0/)
656       
657      ierr = NF_INQ_VARID (nid, "co2ice", nvarid)
658      IF (ierr .NE. NF_NOERR) THEN
659         PRINT*, "lect_start_archive: Le champ <co2ice> est absent"
660         CALL abort
661      ENDIF
662#ifdef NC_DOUBLE
663      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,co2iceold)
664#else
665      ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,co2iceold)
666#endif
667      IF (ierr .NE. NF_NOERR) THEN
668         PRINT*, "lect_start_archive: Lecture echouee pour <co2ice>"
669         PRINT*, NF_STRERROR(ierr)
670         CALL abort
671      ENDIF
672c
673      ierr = NF_INQ_VARID (nid, "emis", nvarid)
674      IF (ierr .NE. NF_NOERR) THEN
675         PRINT*, "lect_start_archive: Le champ <emis> est absent"
676         CALL abort
677      ENDIF
678#ifdef NC_DOUBLE
679      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,emisold)
680#else
681      ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,emisold)
682#endif
683      IF (ierr .NE. NF_NOERR) THEN
684         PRINT*, "lect_start_archive: Lecture echouee pour <emis>"
685         CALL abort
686      ENDIF
687c
688      ierr = NF_INQ_VARID (nid, "ps", nvarid)
689      IF (ierr .NE. NF_NOERR) THEN
690         PRINT*, "lect_start_archive: Le champ <ps> est absent"
691         CALL abort
692      ENDIF
693#ifdef NC_DOUBLE
694      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,psold)
695#else
696      ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,psold)
697#endif
698      IF (ierr .NE. NF_NOERR) THEN
699         PRINT*, "lect_start_archive: Lecture echouee pour <ps>"
700         CALL abort
701      ENDIF
702c
703      ierr = NF_INQ_VARID (nid, "tsurf", nvarid)
704      IF (ierr .NE. NF_NOERR) THEN
705         PRINT*, "lect_start_archive: Le champ <tsurf> est absent"
706         CALL abort
707      ENDIF
708#ifdef NC_DOUBLE
709      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,tsurfold)
710#else
711      ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,tsurfold)
712#endif
713      IF (ierr .NE. NF_NOERR) THEN
714         PRINT*, "lect_start_archive: Lecture echouee pour <tsurf>"
715         CALL abort
716      ENDIF
717c
718      ierr = NF_INQ_VARID (nid, "q2surf", nvarid)
719      IF (ierr .NE. NF_NOERR) THEN
720         PRINT*, "lect_start_archive: Le champ <q2surf> est absent"
721         CALL abort
722      ENDIF
723#ifdef NC_DOUBLE
724      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,q2old)
725#else
726      ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,q2old)
727#endif
728      IF (ierr .NE. NF_NOERR) THEN
729         PRINT*, "lect_start_archive: Lecture echouee pour <q2surf>"
730         CALL abort
731      ENDIF
732c
733      write(*,*)"lect_start_archive: rlonuold:"
734     &           ,rlonuold," rlatvold:",rlatvold
735      write(*,*)
736c
737
738c tracers: the 2 last ones are kept the 2 last one.
739c the others keep their rank. ! No longer true.
740c -------------------------------------------
741! Surface tracers:     
742      do iq=1,nqmx
743        call initial0((jmold+1)*(imold+1), qsurfold(1,1,iq))
744      enddo
745
746      DO iq=1,nqmx
747!      iq=nqold
748!        write(str2,'(i2.2)') iq
749        IF (oldtracernames) THEN
750          txt=" "
751          write(txt,'(a5,i2.2)')'qsurf',iq
752        ELSE
753          txt=trim(tnom(iq))//"_surf"
754          if (txt.eq."h2o_vap") then
755            ! There is no surface tracer for h2o_vap;
756            ! "h2o_ice" should be loaded instead
757            txt="h2o_ice_surf"
758            write(*,*) 'lect_start_archive: loading surface tracer',
759     &                     ' h2o_ice instead of h2o_vap'
760          endif
761        ENDIF ! of IF (oldtracernames)
762        write(*,*) "lect_start_archive: loading tracer ",trim(txt)
763        ierr = NF_INQ_VARID (nid,txt,nvarid)
764        IF (ierr .NE. NF_NOERR) THEN
765          PRINT*, "lect_start_archive: ",
766     &              " Tracer <",trim(txt),"> not found"
767          CALL abort
768        ENDIF
769#ifdef NC_DOUBLE
770        ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,
771     &          qsurfold(1,1,iq))
772#else
773        ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,
774     &          qsurfold(1,1,iq))
775#endif
776        IF (ierr .NE. NF_NOERR) THEN
777          PRINT*, "lect_start_archive: ",
778     &             " Failed loading <",trim(txt),">"
779          write (*,*) trim(txt),'     set to 0'
780          call initial0((jmold+1)*(imold+1), qsurfold(1,1,iq))
781        ENDIF
782
783      ENDDO ! of DO iq=1,nqmx
784
785! The trick below is to read in tracer nqmx-1 (i.e. water ice)
786!      if ((nqold.gt.1).and.(nqmx.gt.1)) then
787!        iq=nqold-1
788!        write(str2,'(i2.2)') iq
789!         ierr = NF_INQ_VARID (nid, "qsurf"//str2, nvarid)
790!         IF (ierr .NE. NF_NOERR) THEN
791!            PRINT*, "lect_start_archive:
792!     .               Le champ <","qsurf"//str2,"> est absent"
793!            CALL abort
794!         ENDIF
795!#ifdef NC_DOUBLE
796!         ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,
797!     .          qsurfold(1,1,nqmx-1))
798!#else
799!         ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,
800!     .          qsurfold(1,1,nqmx-1))
801!#endif
802!         IF (ierr .NE. NF_NOERR) THEN
803!            PRINT*, "lect_start_archive:
804!     .               Lecture echouee pour <","qsurf"//str2,">"
805!            write (*,*) 'qsurf'//str2,' set to 0'
806!            call initial0((jmold+1)*(imold+1), qsurfold(1,1,nqmx-1))
807!         ENDIF
808!      endif
809
810! The trick below is to read in tracers except nqmx & nqmx-1 (water vap.& ice)
811!      if (nqold.gt.2) then
812!       do  iq = 1, nqold-2
813!       if (iq.lt.nqmx-1) then
814!         write(str2,'(i2.2)') iq
815!         ierr = NF_INQ_VARID (nid, "qsurf"//str2, nvarid)
816!         IF (ierr .NE. NF_NOERR) THEN
817!            PRINT*, "lect_start_archive:
818!     .               Le champ <","qsurf"//str2,"> est absent"
819!            CALL abort
820!         ENDIF
821!#ifdef NC_DOUBLE
822!         ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,
823!     .          qsurfold(1,1,iq))
824!#else
825!         ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,
826!     .          qsurfold(1,1,iq))
827!#endif
828!         IF (ierr .NE. NF_NOERR) THEN
829!            PRINT*, "lect_start_archive:
830!     .               Lecture echouee pour <","qsurf"//str2,">"
831!            write (*,*) 'qsurf'//str2,' set to 0'
832!            call initial0((jmold+1)*(imold+1), qsurfold(1,1,iq))
833!         ENDIF
834!       end if
835!       end do
836!      end if
837
838!-----------------------------------------------------------------------
839! 5.2 Read 3D subterranean fields
840!-----------------------------------------------------------------------
841
842      start=(/1,1,1,memo/)
843      count=(/imold+1,jmold+1,nsoilold,1/)
844!
845! Read soil temperatures
846!
847      if (olddepthdef) then ! tsoil stored using the 'old format'
848         start=(/1,1,memo,0/)
849         count=(/imold+1,jmold+1,1,0/) ! because the "Tg" are 2D datasets
850       do isoil=1,nsoilold
851!        write(*,*)'isoil',isoil
852         write(str2,'(i2.2)') isoil
853c
854         ierr = NF_INQ_VARID (nid, "Tg"//str2, nvarid)
855         IF (ierr .NE. NF_NOERR) THEN
856            PRINT*, "lect_start_archive: ",
857     &              "Field <","Tg"//str2,"> not found"
858            CALL abort
859         ENDIF
860#ifdef NC_DOUBLE
861         ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,
862     &          tsoilold(1,1,isoil))
863#else
864         ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,
865     &          tsoilold(1,1,isoil))
866#endif
867         IF (ierr .NE. NF_NOERR) THEN
868            PRINT*, "lect_start_archive: ",
869     &            "Failed reading <","Tg"//str2,">"
870            CALL abort
871         ENDIF
872c
873       enddo ! of do isoil=1,nsoilold
874     
875      ! reset 'start' and 'count' to "3D" behaviour
876      start=(/1,1,1,memo/)
877      count=(/imold+1,jmold+1,nsoilold,1/)
878     
879      else
880       write(*,*) "lect_start_archive: loading tsoil "
881       ierr=NF_INQ_VARID(nid,"tsoil",nvarid)
882       if (ierr.ne.NF_NOERR) then
883        write(*,*)"lect_start_archive: Cannot find <tsoil>"
884        call abort
885       else
886#ifdef NC_DOUBLE
887      ierr=NF_GET_VARA_DOUBLE(nid,nvarid,start,count,tsoilold)
888#else
889      ierr=NF_GET_VARA_REAL(nid,nvarid,start,count,tsoilold)
890#endif
891       endif ! of if (ierr.ne.NF_NOERR)
892       
893      endif ! of if (olddepthdef)
894
895!
896! Read soil thermal inertias
897!
898!      if (.not.olddepthdef) then ! no thermal inertia data in "old" archives
899!       ierr=NF_INQ_VARID(nid,"inertiedat",nvarid)
900!       if (ierr.ne.NF_NOERR) then
901!        write(*,*)"lect_start_archive: Cannot find <inertiedat>"
902!       call abort
903!       else
904!#ifdef NC_DOUBLE
905!      ierr=NF_GET_VARA_DOUBLE(nid,nvarid,start,count,inertiedatold)
906!#else
907!      ierr=NF_GET_VARA_REAL(nid,nvarid,start,count,inertiedatold)
908!#endif
909!       endif ! of if (ierr.ne.NF_NOERR)
910!      endif
911
912c-----------------------------------------------------------------------
913c 5.3   Lecture des champs 3D (t,u,v, q2atm,q)
914c-----------------------------------------------------------------------
915
916      start=(/1,1,1,memo/)
917      count=(/imold+1,jmold+1,lmold,1/)
918
919c
920      ierr = NF_INQ_VARID (nid,"temp", nvarid)
921      IF (ierr .NE. NF_NOERR) THEN
922         PRINT*, "lect_start_archive: Le champ <temp> est absent"
923         CALL abort
924      ENDIF
925#ifdef NC_DOUBLE
926      ierr = NF_GET_VARA_DOUBLE(nid, nvarid, start, count, told)
927#else
928      ierr = NF_GET_VARA_REAL(nid, nvarid, start, count, told)
929#endif
930      IF (ierr .NE. NF_NOERR) THEN
931         PRINT*, "lect_start_archive: Lecture echouee pour <temp>"
932         CALL abort
933      ENDIF
934c
935      ierr = NF_INQ_VARID (nid,"u", nvarid)
936      IF (ierr .NE. NF_NOERR) THEN
937         PRINT*, "lect_start_archive: Le champ <u> est absent"
938         CALL abort
939      ENDIF
940#ifdef NC_DOUBLE
941      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,uold)
942#else
943      ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,uold)
944#endif
945      IF (ierr .NE. NF_NOERR) THEN
946         PRINT*, "lect_start_archive: Lecture echouee pour <u>"
947         CALL abort
948      ENDIF
949c
950      ierr = NF_INQ_VARID (nid,"v", nvarid)
951      IF (ierr .NE. NF_NOERR) THEN
952         PRINT*, "lect_start_archive: Le champ <v> est absent"
953         CALL abort
954      ENDIF
955#ifdef NC_DOUBLE
956      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,vold)
957#else
958      ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,vold)
959#endif
960      IF (ierr .NE. NF_NOERR) THEN
961         PRINT*, "lect_start_archive: Lecture echouee pour <v>"
962         CALL abort
963      ENDIF
964c
965      ierr = NF_INQ_VARID (nid,"q2atm", nvarid)
966      IF (ierr .NE. NF_NOERR) THEN
967         PRINT*, "lect_start_archive: Le champ <q2atm> est absent"
968         CALL abort
969      ENDIF
970#ifdef NC_DOUBLE
971      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,q2old(1,1,2))
972#else
973      ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,q2old(1,1,2))
974#endif
975      IF (ierr .NE. NF_NOERR) THEN
976         PRINT*, "lect_start_archive: Lecture echouee pour <q2atm>"
977         CALL abort
978      ENDIF
979c
980
981c tracers: the 2 last ones are kept the 2 last one.
982c the others keep their rank. ! No longer true.
983c -------------------------------------------
984! Tracers:     
985      do iq=1,nqmx
986         call initial0((jmold+1)*(imold+1)*lmold,qold(1,1,1,iq) )
987      enddo
988
989      DO iq=1,nqmx
990!      iq=nqold
991!        write(str2,'(i2.2)') iq
992        IF (oldtracernames) THEN
993          txt=" "
994          write(txt,'(a1,i2.2)')'q',iq
995        ELSE
996          txt=tnom(iq)
997        ENDIF
998        write(*,*)"lect_start_archive: loading tracer ",trim(txt)
999        ierr = NF_INQ_VARID (nid,txt,nvarid)
1000        IF (ierr .NE. NF_NOERR) THEN
1001            PRINT*, "lect_start_archive: ",
1002     &              " Tracer <",trim(txt),"> not found"
1003            CALL abort
1004        ENDIF
1005#ifdef NC_DOUBLE
1006        ierr=NF_GET_VARA_DOUBLE(nid,nvarid,start,count,qold(1,1,1,iq))
1007#else
1008        ierr=NF_GET_VARA_REAL(nid,nvarid,start,count,qold(1,1,1,iq))
1009#endif
1010        IF (ierr .NE. NF_NOERR) THEN
1011          PRINT*, "lect_start_archive: ",
1012     &             "  Failed loading <",trim(txt),">"
1013          write (*,*) trim(txt),'      set to 1.E-30'
1014          do l=1,lmold
1015            do j=1,jmold+1
1016              do i=1,imold+1
1017                 qold(i,j,l,iq)=1.e-30
1018              end do
1019            end do
1020          end do
1021        ENDIF
1022
1023      ENDDO ! of DO iq=1,nqmx
1024
1025! The trick below is to read in tracer nqmx-1 (i.e. water ice)
1026!      if ((nqold.gt.1).and.(nqmx.gt.1)) then
1027!        iq=nqold-1
1028!        write(str2,'(i2.2)') iq
1029!         ierr = NF_INQ_VARID (nid, "q"//str2, nvarid)
1030!         IF (ierr .NE. NF_NOERR) THEN
1031!            PRINT*, "lect_start_archive:
1032!     .               Le champ <","q"//str2,"> est absent"
1033!            CALL abort
1034!         ENDIF
1035!#ifdef NC_DOUBLE
1036!         ierr= NF_GET_VARA_DOUBLE(nid,nvarid,start,count,
1037!     .                            qold(1,1,1,nqmx-1))
1038!#else
1039!         ierr= NF_GET_VARA_REAL(nid,nvarid,start,count,
1040!     .                            qold(1,1,1,nqmx-1))
1041!#endif
1042!         IF (ierr .NE. NF_NOERR) THEN
1043!            PRINT*, "lect_start_archive:
1044!     .               Lecture echouee pour <","q"//str2,">"
1045!            write (*,*) 'q'//str2,' set to 1.E-30'
1046!            do l=1,lmold
1047!              do j=1,jmold+1
1048!                do i=1,imold+1
1049!                   qold(1,1,1,nqmx-1)=1.e-30
1050!                end do
1051!              end do
1052!            end do
1053!
1054!         ENDIF
1055!      endif
1056
1057! The trick below is to read in tracers except nqmx & nqmx-1 (water vap.& ice)
1058!      if (nqold.gt.2) then
1059!       do  iq = 1, nqold-2
1060!       if (iq.lt.nqmx-1) then
1061!         write(str2,'(i2.2)') iq
1062!         ierr = NF_INQ_VARID (nid, "q"//str2, nvarid)
1063!         IF (ierr .NE. NF_NOERR) THEN
1064!            PRINT*, "lect_start_archive:
1065!     .               Le champ <","q"//str2,"> est absent"
1066!            CALL abort
1067!         ENDIF
1068!#ifdef NC_DOUBLE
1069!         ierr= NF_GET_VARA_DOUBLE(nid,nvarid,start,count,qold(1,1,1,iq))
1070!#else
1071!         ierr= NF_GET_VARA_REAL(nid,nvarid,start,count,qold(1,1,1,iq))
1072!#endif
1073!         IF (ierr .NE. NF_NOERR) THEN
1074!            PRINT*, "lect_start_archive:
1075!     .               Lecture echouee pour <","q"//str2,">"
1076!            write (*,*) 'q'//str2,' set to 1.E-30 '
1077!            do l=1,lmold
1078!              do j=1,jmold+1
1079!                do i=1,imold+1
1080!                   qold(1,1,1,iq)=1.e-30
1081!                end do
1082!              end do
1083!            end do
1084
1085!         ENDIF
1086!       end if
1087!       end do
1088!      end if
1089
1090c Chemin pour trouver les donnees de surface (albedo, relief, th.inertia...)
1091c -------------------------------------------------------------------------
1092
1093!      datapath = '/users/forget/gcm/data_mars_gcm'
1094
1095
1096!=======================================================================
1097! 6. Interpolation from old grid to new grid
1098!=======================================================================
1099
1100c=======================================================================
1101c   INTERPOLATION DANS LA NOUVELLE GRILLE et initialisation des variables
1102c=======================================================================
1103c  Interpolation horizontale puis passage dans la grille physique pour
1104c  les variables physique
1105c  Interpolation verticale puis horizontale pour chaque variable 3D
1106c=======================================================================
1107
1108c-----------------------------------------------------------------------
1109c 6.1   Variable 2d :
1110c-----------------------------------------------------------------------
1111c Relief
1112      call interp_horiz (phisold,phisold_newgrid,imold,jmold,iim,jjm,1,
1113     &                   rlonuold,rlatvold,rlonu,rlatv)
1114
1115c Glace CO2
1116      call interp_horiz (co2iceold,co2ices,imold,jmold,iim,jjm,1,
1117     &                   rlonuold,rlatvold,rlonu,rlatv)
1118
1119c Temperature de surface
1120      call interp_horiz (tsurfold,tsurfs,imold,jmold,iim,jjm,1,
1121     &                   rlonuold,rlatvold,rlonu,rlatv)
1122      call gr_dyn_fi (1,iim+1,jjm+1,ngridmx,tsurfs,tsurf)
1123c     write(44,*) 'tsurf', tsurf
1124
1125c Temperature du sous-sol
1126!      call interp_horiz(tsoilold,tsoils,
1127!     &                  imold,jmold,iim,jjm,nsoilmx,
1128!     &                   rlonuold,rlatvold,rlonu,rlatv)
1129!      call gr_dyn_fi (nsoilmx,iim+1,jjm+1,ngridmx,tsoils,tsoil)
1130c     write(45,*) 'tsoil',tsoil
1131
1132c Emissivite de la surface
1133      call interp_horiz (emisold,emiss,imold,jmold,iim,jjm,1,
1134     &                   rlonuold,rlatvold,rlonu,rlatv)
1135      call gr_dyn_fi (1,iim+1,jjm+1,ngridmx,emiss,emis)
1136c     write(46,*) 'emis',emis
1137c-----------------------------------------------------------------------
1138c 6.1.2 Traitement special de la pression au sol :
1139c-----------------------------------------------------------------------
1140
1141c  Extrapolation la pression dans la nouvelle grille
1142      call interp_horiz(psold,ps,imold,jmold,iim,jjm,1,
1143     &                   rlonuold,rlatvold,rlonu,rlatv)
1144
1145c-----------------------------------------------------------------------
1146c       On assure la conservation de la masse de l'atmosphere + calottes
1147c-----------------------------------------------------------------------
1148
1149      ptotal =  0.
1150      co2icetotal = 0.
1151      DO j=1,jjp1
1152         DO i=1,iim
1153            ptotal=ptotal+ps(i,j)*aire(i,j)/g
1154            co2icetotal = co2icetotal + co2iceS(i,j)*aire(i,j)
1155         END DO
1156      END DO
1157
1158      write(*,*)
1159      write(*,*)'Ancienne grille: masse de l atm :',ptotalold
1160      write(*,*)'Nouvelle grille: masse de l atm :',ptotal
1161      write (*,*) 'Ratio new atm./ old atm =', ptotal/ptotalold
1162      write(*,*)
1163      write(*,*)'Ancienne grille: masse de la glace CO2:',co2icetotalold
1164      write(*,*)'Nouvelle grille: masse de la glace CO2:',co2icetotal
1165      write(*,*)'Ratio new ice./old ice =',co2icetotal/co2icetotalold
1166      write(*,*)
1167
1168
1169      DO j=1,jjp1
1170         DO i=1,iip1
1171            ps(i,j)=ps(i,j) * ptotalold/ptotal
1172         END DO
1173      END DO
1174
1175      if ( co2icetotalold.gt.0.) then
1176         DO j=1,jjp1
1177            DO i=1,iip1
1178               co2iceS(i,j)=co2iceS(i,j) * co2icetotalold/co2icetotal
1179            END DO
1180         END DO
1181      end if
1182
1183c-----------------------------------------------------------------------
1184c 6.2 Subterranean 3d variables:
1185c-----------------------------------------------------------------------
1186
1187c-----------------------------------------------------------------------
1188c 6.2.1 Thermal Inertia
1189c       Note: recall that inertiedat is a common in "comsoil.h"
1190c-----------------------------------------------------------------------
1191
1192      ! depth-wise interpolation, if required
1193      if (depthinterpol.and.(.not.olddepthdef)) then
1194        allocate(oldval(nsoilold))
1195        allocate(newval(nsoilmx))
1196        write(*,*)'lect_start_archive: WARNING: vertical interpolation o
1197     &f soil thermal inertia; might be wiser to reset it.'
1198        write(*,*)
1199       
1200        do i=1,imold+1
1201         do j=1,jmold+1
1202           !copy old values
1203           oldval(1:nsoilold)=inertiedatold(i,j,1:nsoilold)
1204           !interpolate
1205           call interp_line(mlayerold,oldval,nsoilold,
1206     &                     mlayer,newval,nsoilmx)
1207           !copy interpolated values
1208           inertiedatoldnew(i,j,1:nsoilmx)=newval(1:nsoilmx)
1209         enddo
1210        enddo
1211        ! cleanup
1212        deallocate(oldval)
1213        deallocate(newval)
1214      endif !of if (depthinterpol)
1215
1216      if (therminertia_3D) then
1217        ! We have inertiedatold
1218       if((imold.ne.iim).or.(jmold.ne.jjm)) then
1219       write(*,*)'lect_start_archive: WARNING: horizontal interpolation
1220     &of thermal inertia; might be better to reset it.'
1221       write(*,*)
1222       endif
1223       
1224        ! Do horizontal interpolation
1225        if (depthinterpol) then
1226          call interp_horiz(inertiedatoldnew,inertiedatS,
1227     &                  imold,jmold,iim,jjm,nsoilmx,
1228     &                   rlonuold,rlatvold,rlonu,rlatv)
1229        else
1230          call interp_horiz(inertiedatold,inertiedatS,
1231     &                  imold,jmold,iim,jjm,nsoilold,
1232     &                   rlonuold,rlatvold,rlonu,rlatv)
1233        endif ! of if (depthinterpol)
1234
1235      else ! no 3D thermal inertia data
1236       write(*,*)'lect_start_archive: using reference surface inertia'
1237        ! Use surface inertia (and extend it to all depths)
1238        do i=1,nsoilmx
1239         inertiedatS(1:iip1,1:jjp1,i)=surfith(1:iip1,1:jjp1)
1240        enddo
1241        ! Build an old resolution surface thermal inertia
1242        ! (will be needed for tsoil interpolation)
1243        call interp_horiz(surfith,surfithold,
1244     &                    iim,jjm,imold,jmold,1,
1245     &                    rlonu,rlatv,rlonuold,rlatvold)
1246      endif
1247
1248
1249      ! Reshape inertiedatS to scalar grid as inertiedat
1250      call gr_dyn_fi (nsoilmx,iim+1,jjm+1,ngridmx,
1251     &                  inertiedatS,inertiedat)
1252     
1253c-----------------------------------------------------------------------
1254c 6.2.2 Soil temperature
1255c-----------------------------------------------------------------------
1256!      write(*,*) 'Soil'
1257      ! Recast temperatures along soil depth, if necessary
1258      if (olddepthdef) then
1259        allocate(oldgrid(nsoilold+1))
1260        allocate(oldval(nsoilold+1))
1261        allocate(newval(nsoilmx))
1262        do i=1,imold+1
1263         do j=1,jmold+1
1264           ! copy values
1265           oldval(1)=tsurfS(i,j)
1266           oldval(2:nsoilold+1)=tsoilold(i,j,1:nsoilold)
1267           ! build vertical coordinate
1268           oldgrid(1)=0. ! ground
1269           oldgrid(2:nsoilold+1)=mlayerold(1:nsoilold)*
1270     &                (surfithold(i,j)/1.e6)
1271          ! Note; at this stage, we impose volcapa=1.e6 above
1272          ! since volcapa isn't set in old soil definitions
1273
1274          ! interpolate
1275          call interp_line(oldgrid,oldval,nsoilold+1,
1276     &                     mlayer,newval,nsoilmx)
1277         ! copy result in tsoilold
1278         tsoiloldnew(i,j,1:nsoilmx)=newval(1:nsoilmx)
1279         enddo
1280        enddo
1281        ! cleanup
1282        deallocate(oldgrid)
1283        deallocate(oldval)
1284        deallocate(newval)
1285
1286      else
1287       if (depthinterpol) then ! if vertical interpolation is required
1288        allocate(oldgrid(nsoilold+1))
1289        allocate(oldval(nsoilold+1))
1290        allocate(newval(nsoilmx))
1291        ! build vertical coordinate
1292        oldgrid(1)=0. ! ground
1293        oldgrid(2:nsoilold+1)=mlayerold(1:nsoilold)
1294        do i=1,imold+1
1295         do j=1,jmold+1
1296           ! copy values
1297           oldval(1)=tsurfS(i,j)
1298           oldval(2:nsoilold+1)=tsoilold(i,j,1:nsoilold)
1299          ! interpolate
1300          call interp_line(oldgrid,oldval,nsoilold+1,
1301     &                     mlayer,newval,nsoilmx)
1302         ! copy result in tsoilold
1303         tsoiloldnew(i,j,1:nsoilmx)=newval(1:nsoilmx)
1304         enddo
1305        enddo
1306!       write(*,*)'tsoiloldnew(1,1,1):',tsoiloldnew(1,1,1)
1307        ! cleanup
1308        deallocate(oldgrid)
1309        deallocate(oldval)
1310        deallocate(newval)
1311       
1312       else
1313        tsoiloldnew(:,:,:)=tsoilold(:,:,:)
1314       endif ! of if (depthinterpol)
1315      endif ! of if (olddepthdef)
1316
1317      ! Do the horizontal interpolation
1318       call interp_horiz(tsoiloldnew,tsoilS,
1319     &                  imold,jmold,iim,jjm,nsoilmx,
1320     &                   rlonuold,rlatvold,rlonu,rlatv)
1321
1322      ! Reshape tsoilS to scalar grid as tsoil
1323       call gr_dyn_fi (nsoilmx,iim+1,jjm+1,ngridmx,tsoilS,tsoil)
1324
1325
1326
1327c-----------------------------------------------------------------------
1328c 6.3 Variable 3d :
1329c-----------------------------------------------------------------------
1330     
1331c temperatures atmospheriques
1332      write (*,*) 'lect_start_archive: told ', told (1,jmold+1,1)  ! INFO
1333      call interp_vert
1334     &    (told,var,lmold,llm,apsold,bpsold,aps,bps,
1335     &     psold,(imold+1)*(jmold+1))
1336      write (*,*) 'lect_start_archive: var ', var (1,jmold+1,1)  ! INFO
1337      call interp_horiz(var,t,imold,jmold,iim,jjm,llm,
1338     &                   rlonuold,rlatvold,rlonu,rlatv)
1339      write (*,*) 'lect_start_archive: t ', t(1,jjp1,1)  ! INFO
1340
1341c q2 : pbl wind variance
1342      write (*,*) 'lect_start_archive: q2old ', q2old (1,2,1)  ! INFO
1343      call interp_vert (q2old,varp1,lmold+1,llm+1,
1344     &     apsold,bpsold,aps,bps,psold,(imold+1)*(jmold+1))
1345      write (*,*) 'lect_start_archive: varp1 ', varp1 (1,2,1)  ! INFO
1346      call interp_horiz(varp1,q2s,imold,jmold,iim,jjm,llm+1,
1347     &                   rlonuold,rlatvold,rlonu,rlatv)
1348      write (*,*) 'lect_start_archive: q2s ', q2s (1,2,1)  ! INFO
1349      call gr_dyn_fi (llm+1,iim+1,jjm+1,ngridmx,q2s,q2)
1350      write (*,*) 'lect_start_archive: q2 ', q2 (1,2)  ! INFO
1351c     write(47,*) 'q2',q2
1352
1353c calcul des champ de vent; passage en vent covariant
1354      write (*,*) 'lect_start_archive: uold ', uold (1,2,1)  ! INFO
1355      call interp_vert
1356     & (uold,var,lmold,llm,apsold,bpsold,aps,bps,
1357     &  psold,(imold+1)*(jmold+1))
1358      write (*,*) 'lect_start_archive: var ', var (1,2,1)  ! INFO
1359      call interp_horiz(var,us,imold,jmold,iim,jjm,llm,
1360     &                   rlonuold,rlatvold,rlonu,rlatv)
1361      write (*,*) 'lect_start_archive: us ', us (1,2,1)   ! INFO
1362
1363      call interp_vert
1364     & (vold,var,lmold,llm,
1365     &  apsold,bpsold,aps,bps,psold,(imold+1)*(jmold+1))
1366      call interp_horiz(var,vs,imold,jmold,iim,jjm,llm,
1367     &                   rlonuold,rlatvold,rlonu,rlatv)
1368      call scal_wind(us,vs,unat,vnat)
1369      write (*,*) 'lect_start_archive: unat ', unat (1,2,1)    ! INFO
1370      do l=1,llm
1371        do j = 1, jjp1
1372          do i=1,iip1
1373            ucov( i,j,l ) = unat( i,j,l ) * cu(i,j)
1374c           ucov( i,j,l ) = 0
1375          end do
1376        end do
1377      end do
1378      write (*,*) 'lect_start_archive: ucov ', ucov (1,2,1)  ! INFO
1379c     write(48,*) 'ucov',ucov
1380      do l=1,llm
1381        do j = 1, jjm
1382          do i=1,iim
1383            vcov( i,j,l ) = vnat( i,j,l ) * cv(i,j)
1384c           vcov( i,j,l ) = 0
1385          end do
1386          vcov( iip1,j,l ) = vcov( 1,j,l )
1387        end do
1388      end do
1389c     write(49,*) 'ucov',vcov
1390
1391c traceurs surface
1392      do iq = 1, nqmx
1393            call interp_horiz(qsurfold(1,1,iq) ,qsurfs(1,1,iq),
1394     &                  imold,jmold,iim,jjm,1,
1395     &                  rlonuold,rlatvold,rlonu,rlatv)
1396      enddo
1397
1398      call gr_dyn_fi (nqmx,iim+1,jjm+1,ngridmx,qsurfs,qsurf)
1399
1400c traceurs 3D
1401      do  iq = 1, nqmx
1402            call interp_vert(qold(1,1,1,iq),var,lmold,llm,
1403     &        apsold,bpsold,aps,bps,psold,(imold+1)*(jmold+1))
1404            call interp_horiz(var,q(1,1,1,iq),imold,jmold,iim,jjm,llm,
1405     &                  rlonuold,rlatvold,rlonu,rlatv)
1406      enddo
1407cccccccccccccccccccccccccccccc     
1408c  make sure that sum of q = 1     
1409c dominent species is = 1 - sum(all other species)     
1410cccccccccccccccccccccccccccccc     
1411c      iqmax=1
1412c     
1413c      if (nqold.gt.10) then
1414c       do l=1,llm
1415c        do j=1,jjp1
1416c          do i=1,iip1
1417c           do iq=1,nqold
1418c            if (q(i,j,l,iq).gt.q(i,j,l,iqmax)) then
1419c              iqmax=iq
1420c            endif
1421c           enddo
1422c           q(i,j,l,iqmax)=1.
1423c           qtot(i,j,l)=0
1424c           do iq=1,nqold
1425c            if (iq.ne.iqmax) then       
1426c              q(i,j,l,iqmax)=q(i,j,l,iqmax)-q(i,j,l,iq)       
1427c            endif
1428c           enddo !iq
1429c           do iq=1,nqold
1430c            qtot(i,j,l)=qtot(i,j,l)+q(i,j,l,iq)
1431c            if (i.eq.1.and.j.eq.1.and.l.Eq.1) write(*,*)' qtot(i,j,l)',
1432c     $    qtot(i,j,l)
1433c           enddo !iq
1434c          enddo !i   
1435c         enddo !j   
1436c       enddo !l 
1437c      endif
1438ccccccccccccccccccccccccccccccc
1439
1440c     Periodicite :
1441      do  iq = 1, nqmx
1442         do l=1, llm
1443            do j = 1, jjp1
1444               q(iip1,j,l,iq) = q(1,j,l,iq)
1445            end do
1446         end do
1447      enddo
1448     
1449      call gr_dyn_fi (1,iim+1,jjm+1,ngridmx,co2ices,co2ice)
1450
1451c-----------------------------------------------------------------------
1452c   Initialisation  h:  (passage de t -> h)
1453c-----------------------------------------------------------------------
1454
1455      DO l=1,llm
1456         DO j=1,jjp1
1457            DO i=1,iim
1458               h(i,j,l) = t(i,j,l)*((ps(i,j)/preff)**kappa)
1459            END DO
1460            h(iip1,j,l) =  h(1,j,l)
1461         END DO
1462      END DO
1463
1464
1465c***********************************************************************
1466c***********************************************************************
1467c     Fin subroutine lecture ini
1468c***********************************************************************
1469c***********************************************************************
1470
1471      deallocate(timelist)
1472      deallocate(rlonuold, rlatvold)
1473      deallocate(rlonvold, rlatuold)
1474      deallocate(apsold,bpsold)
1475      deallocate(uold)
1476      deallocate(vold)
1477      deallocate(told)
1478      deallocate(psold)
1479      deallocate(phisold)
1480      deallocate(qold)
1481      deallocate(co2iceold)
1482      deallocate(tsurfold)
1483      deallocate(emisold)
1484      deallocate(q2old)
1485      deallocate(tsoilold)
1486      deallocate(tsoiloldnew)
1487      deallocate(inertiedatold)
1488      deallocate(inertiedatoldnew)
1489      deallocate(surfithold)
1490      deallocate(mlayerold)
1491      deallocate(qsurfold)
1492      deallocate(var,varp1)
1493
1494!      write(*,*)'lect_start_archive: END'
1495      return
1496      end
Note: See TracBrowser for help on using the repository browser.