source: trunk/LMDZ.MARS/libf/dyn3d/lect_start_archive.F @ 1036

Last change on this file since 1036 was 1036, checked in by emillour, 12 years ago

Mars GCM: (a first step towards using parallel dynamics)

  • IMPORTANT CHANGE: Implemented dynamic tracers. It is no longer necessary to compile the model with the '-t #' option, number of tracers is simply read from tracer.def file (as before). Adapted makegcm_* scripts (and co.) accordingly. Technical aspects of the switch to dynamic tracers are:
    • advtrac.h (in dyn3d) removed and replaced by module infotrac.F
    • tracer.h (in phymars) removed and replaced by module tracer_mod.F90 (which contains nqmx, the number of tracers, etc. and can be used anywhere in the physics).
  • Included some side cleanups: removed unused files (in dyn3d) anldoppler2.F, anl_mcdstats.F and anl_stats-diag.F, and all the unecessary dimensions.* files in grid/dimension.
  • Checked that changes are clean and that GCM yields identical results (in debug mode) to previous svn version.

EM
File size: 45.0 KB
Line 
1      SUBROUTINE lect_start_archive(nqtot,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      use infotrac, only: tnom
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,intent(in) :: nqtot ! number of advected tracers
52
53c Variables pour les lectures des fichiers "ini"
54c--------------------------------------------------
55!      INTEGER sizei,
56      integer timelen,dimid
57!      INTEGER length
58!      parameter (length = 100)
59      INTEGER tab0
60      INTEGER isoil,iq,iqmax
61      CHARACTER*2   str2
62
63!      REAL dimfirst(4) ! tableau contenant les 1ers elements des dimensions
64
65!      REAL dimlast(4) ! tableau contenant les derniers elements des dimensions
66
67!      REAL dimcycl(4) ! tableau contenant les periodes des dimensions
68!      CHARACTER*120 dimsource
69!      CHARACTER*16 dimname
70!      CHARACTER*80 dimtitle
71!      CHARACTER*40 dimunits
72!      INTEGER   dimtype
73
74!      INTEGER dimord(4)  ! tableau contenant l''ordre
75!      data dimord /1,2,3,4/ ! de sortie des dimensions
76
77!      INTEGER vardim(4)
78      REAL date
79      INTEGER   memo
80!      character (len=50) :: tmpname
81
82c Variable histoire
83c------------------
84      REAL vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm) ! vents covariants
85      REAL h(iip1,jjp1,llm),ps(iip1,jjp1)
86      REAL q(iip1,jjp1,llm,nqtot),qtot(iip1,jjp1,llm)
87
88c autre variables dynamique nouvelle grille
89c------------------------------------------
90
91c!-*-
92!      integer klatdat,klongdat
93!      PARAMETER (klatdat=180,klongdat=360)
94
95c Physique sur grille scalaire
96c----------------------------
97
98c variable physique
99c------------------
100      REAL tsurf(ngridmx) ! surface temperature
101      REAL tsoil(ngridmx,nsoilmx) ! soil temperature
102      REAL co2ice(ngridmx) ! CO2 ice layer
103      REAL emis(ngridmx)
104      REAL q2(ngridmx,nlayermx+1),qsurf(ngridmx,nqtot)
105c     REAL phisfi(ngridmx)
106
107      INTEGER i,j,l
108      INTEGER nid,nvarid
109c     REAL year_day,periheli,aphelie,peri_day
110c     REAL obliquit,z0,emin_turb,lmixmin
111c     REAL emissiv,emisice(2),albedice(2),tauvis
112c     REAL iceradius(2) , dtemisice(2)
113
114!      EXTERNAL RAN1
115!      REAL RAN1
116!      EXTERNAL geopot,inigeom
117      integer ierr
118!      integer ismin
119!      external ismin
120!      CHARACTER*80 datapath
121      integer, dimension(4) :: start,count
122
123c Variable nouvelle grille naturelle au point scalaire
124c------------------------------------------------------
125      real us(iip1,jjp1,llm),vs(iip1,jjp1,llm)
126      REAL phisold_newgrid(iip1,jjp1)
127      REAL t(iip1,jjp1,llm)
128      real tsurfS(iip1,jjp1),tsoilS(iip1,jjp1,nsoilmx)
129      real inertiedatS(iip1,jjp1,nsoilmx)
130      real co2iceS(iip1,jjp1),emisS(iip1,jjp1)
131      REAL q2S(iip1,jjp1,llm+1),qsurfS(iip1,jjp1,nqtot)
132
133      real ptotal, co2icetotal
134
135c Var intermediaires : vent naturel, mais pas coord scalaire
136c-----------------------------------------------------------
137      real vnat(iip1,jjm,llm),unat(iip1,jjp1,llm)
138
139
140c Variable de l'ancienne grille
141c---------------------------------------------------------
142
143      real, dimension(:), allocatable :: timelist
144      real, dimension(:), allocatable :: rlonuold, rlatvold
145      real, dimension(:), allocatable :: rlonvold, rlatuold
146      real, dimension(:), allocatable :: apsold,bpsold
147      real, dimension(:), allocatable :: mlayerold
148      real, dimension(:,:,:), allocatable :: uold,vold,told,q2old
149      real, dimension(:,:,:), allocatable :: tsoilold,qsurfold
150      real, dimension(:,:,:),allocatable :: tsoiloldnew
151! tsoiloldnew: old soil values, but along new subterranean grid
152      real, dimension(:,:,:), allocatable :: inertiedatold
153! inertiedatoldnew: old inertia values, but along new subterranean grid
154      real, dimension(:,:,:), allocatable :: inertiedatoldnew
155      real, dimension(:,:), allocatable :: psold,phisold
156      real, dimension(:,:), allocatable :: co2iceold,tsurfold
157      real, dimension(:,:), allocatable :: emisold
158      real, dimension(:,:,:,:), allocatable :: qold
159
160      real tab_cntrl(100)
161
162      real ptotalold, co2icetotalold
163
164      logical :: olddepthdef=.false. ! flag
165! olddepthdef=.true. if soil depths are in 'old' (unspecified) format
166      logical :: depthinterpol=.false. ! flag
167! depthinterpol=.true. if interpolation will be requiered
168      logical :: therminertia_3D=.true. ! flag
169! therminertia_3D=.true. if thermal inertia is 3D and read from datafile
170c Variable intermediaires iutilise pour l'extrapolation verticale
171c----------------------------------------------------------------
172      real, dimension(:,:,:), allocatable :: var,varp1
173      real, dimension(:), allocatable :: oldgrid, oldval
174      real, dimension(:), allocatable :: newval
175!      real, dimension(:), allocatable :: oldmlayer
176
177      real surfith(iip1,jjp1) ! surface thermal inertia
178!      real surfithfi(ngridmx)
179      ! surface thermal inertia at old horizontal grid resolution
180      real, dimension(:,:), allocatable :: surfithold
181
182! flag which identifies if archive file is using old tracer names (qsurf01,...)
183      logical :: oldtracernames=.false.
184      integer :: counter
185      character(len=30) :: txt ! to store some text
186      real :: tmpval ! to store a temporary variable/value
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,nqtot
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.nqtot) 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,nqtot))
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,nqtot))
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: <controle> is missing"
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: Failed loading <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: <rlonv> is missing"
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: Failed loading <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: <rlatu> is missing"
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: Failed loading <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: <rlonu> is missing"
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: Failed loading <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: <rlatv> is missing"
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: Failed loading <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: <aps> is missing"
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: Failed loading <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: <bps> is missing"
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: Failed loading <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       ! Also build (default) soil interlayer depth layer(:)
525       do isoil=1,nsoilmx
526         layer(isoil)=sqrt(mlayer(0)*mlayer(1))*
527     &                      ((mlayer(1)/mlayer(0))**(isoil-1))
528       enddo
529      endif ! of if (depthinterpol)
530
531c-----------------------------------------------------------------------
532c 3.5 Read Soil thermal inertia
533c-----------------------------------------------------------------------
534
535      ierr=NF_INQ_VARID(nid,"inertiedat",nvarid)
536      if (ierr.ne.NF_NOERR) then
537       write(*,*)'lect_start_archive: Could not find <inertiedat>'
538       write(*,*)' => Assuming this is an archive in old format'
539       therminertia_3D=.false.
540       write(*,*)' => Thermal inertia will be read from reference file'
541       volcapa=1.e6
542       write(*,*)'    and soil volumetric heat capacity is set to ',
543     &           volcapa
544      else
545#ifdef NC_DOUBLE
546        ierr = NF_GET_VAR_DOUBLE(nid,nvarid,inertiedatold)
547#else
548        ierr = NF_GET_VAR_REAL(nid,nvarid,inertiedatold)
549#endif
550       if (ierr .NE. NF_NOERR) then
551         PRINT*, "lect_start_archive: Failed reading <inertiedat>"
552         CALL abort
553       endif
554      endif
555
556c-----------------------------------------------------------------------
557c 3.6 Lecture geopotentiel au sol
558c-----------------------------------------------------------------------
559c
560      ierr = NF_INQ_VARID (nid, "phisinit", nvarid)
561      IF (ierr .NE. NF_NOERR) THEN
562         PRINT*, "lect_start_archive: <phisinit> is missing"
563         CALL abort
564      ENDIF
565#ifdef NC_DOUBLE
566      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, phisold)
567#else
568      ierr = NF_GET_VAR_REAL(nid, nvarid, phisold)
569#endif
570      IF (ierr .NE. NF_NOERR) THEN
571         PRINT*, "lect_start_archive: Failed loading <phisinit>"
572         CALL abort
573      ENDIF
574
575C-----------------------------------------------------------------------
576c   lecture de "ptotalold" et "co2icetotalold"
577c-----------------------------------------------------------------------
578      ptotalold = tab_cntrl(tab0+49)
579      co2icetotalold = tab_cntrl(tab0+50)
580 
581c-----------------------------------------------------------------------
582c 4. Lecture du temps et choix
583c-----------------------------------------------------------------------
584 
585c  lecture du temps
586c
587      ierr = NF_INQ_DIMID (nid, "Time", nvarid)
588      IF (ierr .NE. NF_NOERR) THEN
589         ierr = NF_INQ_DIMID (nid, "temps", nvarid)
590         IF (ierr .NE. NF_NOERR) THEN
591            PRINT*, "lect_start_archive: <Time> is missing"
592            CALL abort
593         endif
594      ENDIF
595
596      ierr = NF_INQ_VARID (nid, "Time", nvarid)
597      IF (ierr .NE. NF_NOERR) THEN
598         ierr = NF_INQ_VARID (nid, "temps", nvarid)
599      endif
600#ifdef NC_DOUBLE
601      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, timelist)
602#else
603      ierr = NF_GET_VAR_REAL(nid, nvarid, timelist)
604#endif
605      IF (ierr .NE. NF_NOERR) THEN
606         PRINT*, "lect_start_archive: Failed loading <Time>"
607         CALL abort
608      ENDIF
609c
610      write(*,*)
611      write(*,*)
612      write(*,*) 'Dates of the stored initial states:'
613      write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
614      pi=2.*ASIN(1.)
615      do i=1,timelen
616c       call solarlong(timelist(i),sollong(i))
617c       sollong(i) = sollong(i)*180./pi
618c        write(*,*) 'initial state at martian day: ',int(timelist(i))
619        write(*,*) 'initial state at martian day: ',timelist(i)
620c       write(*,6) nint(timelist(i)),nint(mod(timelist(i),669)),
621c    .    sollong(i)
622      end do
623
624   6  FORMAT(i7,i7,f9.3)
625 
626      write(*,*)
627      write(*,*) 'Choose the martian day to use'
628 123  read(*,*,iostat=ierr) date
629      if(ierr.ne.0) goto 123
630      memo = 0
631      do i=1,timelen
632c        if (date.eq.int(timelist(i))) then
633        if (abs(date-timelist(i)).lt.0.01) then
634            memo = i
635        endif
636      end do
637 
638      if (memo.eq.0) then
639        write(*,*)
640        write(*,*)
641        write(*,*) 'Wrong value for day number !!'
642        write(*,*)
643        write(*,*) 'Dates of the stored initial states:'
644        write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
645        do i=1,timelen
646          write(*,*) 'initial state at martian day: ',timelist(i)
647c         write(*,6) nint(timelist(i)),nint(mod(timelist(i),669))
648        end do
649        goto 123
650      endif
651     
652!-----------------------------------------------------------------------
653! 5. Read (time-dependent) data from datafile
654!-----------------------------------------------------------------------
655
656
657c-----------------------------------------------------------------------
658c 5.1 Lecture des champs 2D (co2ice, emis,ps,tsurf,Tg[10], q2surf)
659c-----------------------------------------------------------------------
660 
661      start=(/1,1,memo,0/)
662      count=(/imold+1,jmold+1,1,0/)
663       
664      ierr = NF_INQ_VARID (nid, "co2ice", nvarid)
665      IF (ierr .NE. NF_NOERR) THEN
666         PRINT*, "lect_start_archive: <co2ice> is missing"
667         CALL abort
668      ENDIF
669#ifdef NC_DOUBLE
670      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,co2iceold)
671#else
672      ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,co2iceold)
673#endif
674      IF (ierr .NE. NF_NOERR) THEN
675         PRINT*, "lect_start_archive: Failed loading <co2ice>"
676         PRINT*, NF_STRERROR(ierr)
677         CALL abort
678      ENDIF
679c
680      ierr = NF_INQ_VARID (nid, "emis", nvarid)
681      IF (ierr .NE. NF_NOERR) THEN
682         PRINT*, "lect_start_archive: <emis> is missing"
683         CALL abort
684      ENDIF
685#ifdef NC_DOUBLE
686      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,emisold)
687#else
688      ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,emisold)
689#endif
690      IF (ierr .NE. NF_NOERR) THEN
691         PRINT*, "lect_start_archive: Failed loading <emis>"
692         CALL abort
693      ENDIF
694c
695      ierr = NF_INQ_VARID (nid, "ps", nvarid)
696      IF (ierr .NE. NF_NOERR) THEN
697         PRINT*, "lect_start_archive: <ps> is missing"
698         CALL abort
699      ENDIF
700#ifdef NC_DOUBLE
701      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,psold)
702#else
703      ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,psold)
704#endif
705      IF (ierr .NE. NF_NOERR) THEN
706         PRINT*, "lect_start_archive: Failed loading <ps>"
707         CALL abort
708      ENDIF
709c
710      ierr = NF_INQ_VARID (nid, "tsurf", nvarid)
711      IF (ierr .NE. NF_NOERR) THEN
712         PRINT*, "lect_start_archive: <tsurf> is missing"
713         CALL abort
714      ENDIF
715#ifdef NC_DOUBLE
716      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,tsurfold)
717#else
718      ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,tsurfold)
719#endif
720      IF (ierr .NE. NF_NOERR) THEN
721         PRINT*, "lect_start_archive: Failed loading <tsurf>"
722         CALL abort
723      ENDIF
724c
725      ierr = NF_INQ_VARID (nid, "q2surf", nvarid)
726      IF (ierr .NE. NF_NOERR) THEN
727         PRINT*, "lect_start_archive: <q2surf> is missing"
728         CALL abort
729      ENDIF
730#ifdef NC_DOUBLE
731      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,q2old)
732#else
733      ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,q2old)
734#endif
735      IF (ierr .NE. NF_NOERR) THEN
736         PRINT*, "lect_start_archive: Failed loading <q2surf>"
737         CALL abort
738      ENDIF
739c
740      write(*,*)"lect_start_archive: rlonuold:"
741     &           ,rlonuold," rlatvold:",rlatvold
742      write(*,*)
743c
744
745c tracers: the 2 last ones are kept the 2 last one.
746c the others keep their rank. ! No longer true.
747c -------------------------------------------
748! Surface tracers:     
749      qsurfold(1:imold+1,1:jmold+1,1:nqtot)=0
750
751      DO iq=1,nqtot
752        IF (oldtracernames) THEN
753          txt=" "
754          write(txt,'(a5,i2.2)')'qsurf',iq
755        ELSE
756          txt=trim(tnom(iq))//"_surf"
757          if (txt.eq."h2o_vap") then
758            ! There is no surface tracer for h2o_vap;
759            ! "h2o_ice" should be loaded instead
760            txt="h2o_ice_surf"
761            write(*,*) 'lect_start_archive: loading surface tracer',
762     &                     ' h2o_ice instead of h2o_vap'
763          endif
764        ENDIF ! of IF (oldtracernames)
765        write(*,*) "lect_start_archive: loading tracer ",trim(txt)
766        ierr = NF_INQ_VARID (nid,txt,nvarid)
767        IF (ierr .NE. NF_NOERR) THEN
768          PRINT*, "lect_start_archive: ",
769     &              " Tracer <",trim(txt),"> not found"
770          print*, "which (constant) value should it be initialized to?"
771          read(*,*) tmpval
772          qsurfold(1:imold+1,1:jmold+1,iq)=tmpval
773        ELSE ! tracer exists in file, load it
774#ifdef NC_DOUBLE
775          ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,
776     &          qsurfold(1,1,iq))
777#else
778          ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,
779     &          qsurfold(1,1,iq))
780#endif
781          IF (ierr .NE. NF_NOERR) THEN
782            PRINT*, "lect_start_archive: ",
783     &             " Failed loading <",trim(txt),">"
784            stop
785          ENDIF
786        ENDIF
787
788      ENDDO ! of DO iq=1,nqtot
789
790!-----------------------------------------------------------------------
791! 5.2 Read 3D subterranean fields
792!-----------------------------------------------------------------------
793
794      start=(/1,1,1,memo/)
795      count=(/imold+1,jmold+1,nsoilold,1/)
796!
797! Read soil temperatures
798!
799      if (olddepthdef) then ! tsoil stored using the 'old format'
800         start=(/1,1,memo,0/)
801         count=(/imold+1,jmold+1,1,0/) ! because the "Tg" are 2D datasets
802       do isoil=1,nsoilold
803!        write(*,*)'isoil',isoil
804         write(str2,'(i2.2)') isoil
805c
806         ierr = NF_INQ_VARID (nid, "Tg"//str2, nvarid)
807         IF (ierr .NE. NF_NOERR) THEN
808            PRINT*, "lect_start_archive: ",
809     &              "Field <","Tg"//str2,"> not found"
810            CALL abort
811         ENDIF
812#ifdef NC_DOUBLE
813         ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,
814     &          tsoilold(1,1,isoil))
815#else
816         ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,
817     &          tsoilold(1,1,isoil))
818#endif
819         IF (ierr .NE. NF_NOERR) THEN
820            PRINT*, "lect_start_archive: ",
821     &            "Failed reading <","Tg"//str2,">"
822            CALL abort
823         ENDIF
824c
825       enddo ! of do isoil=1,nsoilold
826     
827      ! reset 'start' and 'count' to "3D" behaviour
828      start=(/1,1,1,memo/)
829      count=(/imold+1,jmold+1,nsoilold,1/)
830     
831      else
832       write(*,*) "lect_start_archive: loading tsoil "
833       ierr=NF_INQ_VARID(nid,"tsoil",nvarid)
834       if (ierr.ne.NF_NOERR) then
835        write(*,*)"lect_start_archive: Cannot find <tsoil>"
836        call abort
837       else
838#ifdef NC_DOUBLE
839      ierr=NF_GET_VARA_DOUBLE(nid,nvarid,start,count,tsoilold)
840#else
841      ierr=NF_GET_VARA_REAL(nid,nvarid,start,count,tsoilold)
842#endif
843       endif ! of if (ierr.ne.NF_NOERR)
844       
845      endif ! of if (olddepthdef)
846
847!
848! Read soil thermal inertias
849!
850!      if (.not.olddepthdef) then ! no thermal inertia data in "old" archives
851!       ierr=NF_INQ_VARID(nid,"inertiedat",nvarid)
852!       if (ierr.ne.NF_NOERR) then
853!        write(*,*)"lect_start_archive: Cannot find <inertiedat>"
854!       call abort
855!       else
856!#ifdef NC_DOUBLE
857!      ierr=NF_GET_VARA_DOUBLE(nid,nvarid,start,count,inertiedatold)
858!#else
859!      ierr=NF_GET_VARA_REAL(nid,nvarid,start,count,inertiedatold)
860!#endif
861!       endif ! of if (ierr.ne.NF_NOERR)
862!      endif
863
864c-----------------------------------------------------------------------
865c 5.3   Lecture des champs 3D (t,u,v, q2atm,q)
866c-----------------------------------------------------------------------
867
868      start=(/1,1,1,memo/)
869      count=(/imold+1,jmold+1,lmold,1/)
870
871c
872      ierr = NF_INQ_VARID (nid,"temp", nvarid)
873      IF (ierr .NE. NF_NOERR) THEN
874         PRINT*, "lect_start_archive: <temp> is missing"
875         CALL abort
876      ENDIF
877#ifdef NC_DOUBLE
878      ierr = NF_GET_VARA_DOUBLE(nid, nvarid, start, count, told)
879#else
880      ierr = NF_GET_VARA_REAL(nid, nvarid, start, count, told)
881#endif
882      IF (ierr .NE. NF_NOERR) THEN
883         PRINT*, "lect_start_archive: Failed loading <temp>"
884         CALL abort
885      ENDIF
886c
887      ierr = NF_INQ_VARID (nid,"u", nvarid)
888      IF (ierr .NE. NF_NOERR) THEN
889         PRINT*, "lect_start_archive: <u> is missing"
890         CALL abort
891      ENDIF
892#ifdef NC_DOUBLE
893      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,uold)
894#else
895      ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,uold)
896#endif
897      IF (ierr .NE. NF_NOERR) THEN
898         PRINT*, "lect_start_archive: Failed loading <u>"
899         CALL abort
900      ENDIF
901c
902      ierr = NF_INQ_VARID (nid,"v", nvarid)
903      IF (ierr .NE. NF_NOERR) THEN
904         PRINT*, "lect_start_archive: <v> is missing"
905         CALL abort
906      ENDIF
907#ifdef NC_DOUBLE
908      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,vold)
909#else
910      ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,vold)
911#endif
912      IF (ierr .NE. NF_NOERR) THEN
913         PRINT*, "lect_start_archive: Failed loading <v>"
914         CALL abort
915      ENDIF
916c
917      ierr = NF_INQ_VARID (nid,"q2atm", nvarid)
918      IF (ierr .NE. NF_NOERR) THEN
919         PRINT*, "lect_start_archive: <q2atm> is missing"
920         CALL abort
921      ENDIF
922#ifdef NC_DOUBLE
923      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,q2old(1,1,2))
924#else
925      ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,q2old(1,1,2))
926#endif
927      IF (ierr .NE. NF_NOERR) THEN
928         PRINT*, "lect_start_archive: Failed loading <q2atm>"
929         CALL abort
930      ENDIF
931c
932
933c tracers: the 2 last ones are kept the 2 last one.
934c the others keep their rank. ! No longer true.
935c -------------------------------------------
936! Tracers:
937      qold(1:imold+1,1:jmold+1,1:lmold,1:nqtot)=0
938
939      DO iq=1,nqtot
940        IF (oldtracernames) THEN
941          txt=" "
942          write(txt,'(a1,i2.2)')'q',iq
943        ELSE
944          txt=tnom(iq)
945        ENDIF
946        write(*,*)"lect_start_archive: loading tracer ",trim(txt)
947        ierr = NF_INQ_VARID (nid,txt,nvarid)
948        IF (ierr .NE. NF_NOERR) THEN
949            PRINT*, "lect_start_archive: ",
950     &              " Tracer <",trim(txt),"> not found"
951          print*, "which (constant) value should it be initialized to?"
952          read(*,*) tmpval
953          qold(1:imold+1,1:jmold+1,1:lmold,iq)=tmpval
954        ELSE ! tracer exists in file, load it
955#ifdef NC_DOUBLE
956         ierr=NF_GET_VARA_DOUBLE(nid,nvarid,start,count,qold(1,1,1,iq))
957#else
958         ierr=NF_GET_VARA_REAL(nid,nvarid,start,count,qold(1,1,1,iq))
959#endif
960          IF (ierr .NE. NF_NOERR) THEN
961            PRINT*, "lect_start_archive: ",
962     &             "  Failed loading <",trim(txt),">"
963            stop
964          ENDIF
965        ENDIF
966
967      ENDDO ! of DO iq=1,nqtot
968
969
970c Chemin pour trouver les donnees de surface (albedo, relief, th.inertia...)
971c -------------------------------------------------------------------------
972
973!      datapath = '/users/forget/gcm/data_mars_gcm'
974
975
976!=======================================================================
977! 6. Interpolation from old grid to new grid
978!=======================================================================
979
980c=======================================================================
981c   INTERPOLATION DANS LA NOUVELLE GRILLE et initialisation des variables
982c=======================================================================
983c  Interpolation horizontale puis passage dans la grille physique pour
984c  les variables physique
985c  Interpolation verticale puis horizontale pour chaque variable 3D
986c=======================================================================
987
988c-----------------------------------------------------------------------
989c 6.1   Variable 2d :
990c-----------------------------------------------------------------------
991c Relief
992      call interp_horiz (phisold,phisold_newgrid,imold,jmold,iim,jjm,1,
993     &                   rlonuold,rlatvold,rlonu,rlatv)
994
995c Glace CO2
996      call interp_horiz (co2iceold,co2ices,imold,jmold,iim,jjm,1,
997     &                   rlonuold,rlatvold,rlonu,rlatv)
998
999c Temperature de surface
1000      call interp_horiz (tsurfold,tsurfs,imold,jmold,iim,jjm,1,
1001     &                   rlonuold,rlatvold,rlonu,rlatv)
1002      call gr_dyn_fi (1,iim+1,jjm+1,ngridmx,tsurfs,tsurf)
1003c     write(44,*) 'tsurf', tsurf
1004
1005c Temperature du sous-sol
1006!      call interp_horiz(tsoilold,tsoils,
1007!     &                  imold,jmold,iim,jjm,nsoilmx,
1008!     &                   rlonuold,rlatvold,rlonu,rlatv)
1009!      call gr_dyn_fi (nsoilmx,iim+1,jjm+1,ngridmx,tsoils,tsoil)
1010c     write(45,*) 'tsoil',tsoil
1011
1012c Emissivite de la surface
1013      call interp_horiz (emisold,emiss,imold,jmold,iim,jjm,1,
1014     &                   rlonuold,rlatvold,rlonu,rlatv)
1015      call gr_dyn_fi (1,iim+1,jjm+1,ngridmx,emiss,emis)
1016c     write(46,*) 'emis',emis
1017c-----------------------------------------------------------------------
1018c 6.1.2 Traitement special de la pression au sol :
1019c-----------------------------------------------------------------------
1020
1021c  Extrapolation la pression dans la nouvelle grille
1022      call interp_horiz(psold,ps,imold,jmold,iim,jjm,1,
1023     &                   rlonuold,rlatvold,rlonu,rlatv)
1024
1025c-----------------------------------------------------------------------
1026c       On assure la conservation de la masse de l'atmosphere + calottes
1027c-----------------------------------------------------------------------
1028
1029      ptotal =  0.
1030      co2icetotal = 0.
1031      DO j=1,jjp1
1032         DO i=1,iim
1033            ptotal=ptotal+ps(i,j)*aire(i,j)/g
1034            co2icetotal = co2icetotal + co2iceS(i,j)*aire(i,j)
1035         END DO
1036      END DO
1037
1038      write(*,*)
1039      write(*,*)'Old grid: mass of the atmosphere :',ptotalold
1040      write(*,*)'New grid: mass of the atmosphere :',ptotal
1041      write (*,*) 'Ratio new atm / old atm =', ptotal/ptotalold
1042      write(*,*)
1043      write(*,*)'Old grid: mass of CO2 ice:',co2icetotalold
1044      write(*,*)'New grid: mass of CO2 ice:',co2icetotal
1045      if (co2icetotalold.ne.0.) then
1046       write(*,*)'Ratio new ice / old ice =',co2icetotal/co2icetotalold
1047      endif
1048      write(*,*)
1049
1050
1051      DO j=1,jjp1
1052         DO i=1,iip1
1053            ps(i,j)=ps(i,j) * ptotalold/ptotal
1054         END DO
1055      END DO
1056
1057      if ( co2icetotalold.gt.0.) then
1058         DO j=1,jjp1
1059            DO i=1,iip1
1060               co2iceS(i,j)=co2iceS(i,j) * co2icetotalold/co2icetotal
1061            END DO
1062         END DO
1063      end if
1064
1065c-----------------------------------------------------------------------
1066c 6.2 Subterranean 3d variables:
1067c-----------------------------------------------------------------------
1068
1069c-----------------------------------------------------------------------
1070c 6.2.1 Thermal Inertia
1071c       Note: recall that inertiedat is a common in "comsoil.h"
1072c-----------------------------------------------------------------------
1073
1074      ! depth-wise interpolation, if required
1075      if (depthinterpol.and.(.not.olddepthdef)) then
1076        allocate(oldval(nsoilold))
1077        allocate(newval(nsoilmx))
1078        write(*,*)'lect_start_archive: WARNING: vertical interpolation o
1079     &f soil thermal inertia; might be wiser to reset it.'
1080        write(*,*)
1081       
1082        do i=1,imold+1
1083         do j=1,jmold+1
1084           !copy old values
1085           oldval(1:nsoilold)=inertiedatold(i,j,1:nsoilold)
1086           !interpolate
1087           call interp_line(mlayerold,oldval,nsoilold,
1088     &                     mlayer,newval,nsoilmx)
1089           !copy interpolated values
1090           inertiedatoldnew(i,j,1:nsoilmx)=newval(1:nsoilmx)
1091         enddo
1092        enddo
1093        ! cleanup
1094        deallocate(oldval)
1095        deallocate(newval)
1096      endif !of if (depthinterpol)
1097
1098      if (therminertia_3D) then
1099        ! We have inertiedatold
1100       if((imold.ne.iim).or.(jmold.ne.jjm)) then
1101       write(*,*)'lect_start_archive: WARNING: horizontal interpolation
1102     &of thermal inertia; might be better to reset it.'
1103       write(*,*)
1104       endif
1105       
1106        ! Do horizontal interpolation
1107        if (depthinterpol) then
1108          call interp_horiz(inertiedatoldnew,inertiedatS,
1109     &                  imold,jmold,iim,jjm,nsoilmx,
1110     &                   rlonuold,rlatvold,rlonu,rlatv)
1111        else
1112          call interp_horiz(inertiedatold,inertiedatS,
1113     &                  imold,jmold,iim,jjm,nsoilold,
1114     &                   rlonuold,rlatvold,rlonu,rlatv)
1115        endif ! of if (depthinterpol)
1116
1117      else ! no 3D thermal inertia data
1118       write(*,*)'lect_start_archive: using reference surface inertia'
1119        ! Use surface inertia (and extend it to all depths)
1120        do i=1,nsoilmx
1121         inertiedatS(1:iip1,1:jjp1,i)=surfith(1:iip1,1:jjp1)
1122        enddo
1123        ! Build an old resolution surface thermal inertia
1124        ! (will be needed for tsoil interpolation)
1125        call interp_horiz(surfith,surfithold,
1126     &                    iim,jjm,imold,jmold,1,
1127     &                    rlonu,rlatv,rlonuold,rlatvold)
1128      endif
1129
1130
1131      ! Reshape inertiedatS to scalar grid as inertiedat
1132      call gr_dyn_fi (nsoilmx,iim+1,jjm+1,ngridmx,
1133     &                  inertiedatS,inertiedat)
1134     
1135c-----------------------------------------------------------------------
1136c 6.2.2 Soil temperature
1137c-----------------------------------------------------------------------
1138!      write(*,*) 'Soil'
1139      ! Recast temperatures along soil depth, if necessary
1140      if (olddepthdef) then
1141        allocate(oldgrid(nsoilold+1))
1142        allocate(oldval(nsoilold+1))
1143        allocate(newval(nsoilmx))
1144        do i=1,imold+1
1145         do j=1,jmold+1
1146           ! copy values
1147           oldval(1)=tsurfold(i,j)
1148           oldval(2:nsoilold+1)=tsoilold(i,j,1:nsoilold)
1149           ! build vertical coordinate
1150           oldgrid(1)=0. ! ground
1151           oldgrid(2:nsoilold+1)=mlayerold(1:nsoilold)*
1152     &                (surfithold(i,j)/1.e6)
1153          ! Note; at this stage, we impose volcapa=1.e6 above
1154          ! since volcapa isn't set in old soil definitions
1155
1156          ! interpolate
1157          call interp_line(oldgrid,oldval,nsoilold+1,
1158     &                     mlayer,newval,nsoilmx)
1159         ! copy result in tsoilold
1160         tsoiloldnew(i,j,1:nsoilmx)=newval(1:nsoilmx)
1161         enddo
1162        enddo
1163        ! cleanup
1164        deallocate(oldgrid)
1165        deallocate(oldval)
1166        deallocate(newval)
1167
1168      else
1169       if (depthinterpol) then ! if vertical interpolation is required
1170        allocate(oldgrid(nsoilold+1))
1171        allocate(oldval(nsoilold+1))
1172        allocate(newval(nsoilmx))
1173        ! build vertical coordinate
1174        oldgrid(1)=0. ! ground
1175        oldgrid(2:nsoilold+1)=mlayerold(1:nsoilold)
1176        do i=1,imold+1
1177         do j=1,jmold+1
1178           ! copy values
1179           oldval(1)=tsurfold(i,j)
1180           oldval(2:nsoilold+1)=tsoilold(i,j,1:nsoilold)
1181          ! interpolate
1182          call interp_line(oldgrid,oldval,nsoilold+1,
1183     &                     mlayer,newval,nsoilmx)
1184         ! copy result in tsoilold
1185         tsoiloldnew(i,j,1:nsoilmx)=newval(1:nsoilmx)
1186         enddo
1187        enddo
1188!       write(*,*)'tsoiloldnew(1,1,1):',tsoiloldnew(1,1,1)
1189        ! cleanup
1190        deallocate(oldgrid)
1191        deallocate(oldval)
1192        deallocate(newval)
1193       
1194       else
1195        tsoiloldnew(:,:,:)=tsoilold(:,:,:)
1196       endif ! of if (depthinterpol)
1197      endif ! of if (olddepthdef)
1198
1199      ! Do the horizontal interpolation
1200       call interp_horiz(tsoiloldnew,tsoilS,
1201     &                  imold,jmold,iim,jjm,nsoilmx,
1202     &                   rlonuold,rlatvold,rlonu,rlatv)
1203
1204      ! Reshape tsoilS to scalar grid as tsoil
1205       call gr_dyn_fi (nsoilmx,iim+1,jjm+1,ngridmx,tsoilS,tsoil)
1206
1207
1208
1209c-----------------------------------------------------------------------
1210c 6.3 Variable 3d :
1211c-----------------------------------------------------------------------
1212     
1213c temperatures atmospheriques
1214      write (*,*) 'lect_start_archive: told ', told (1,jmold+1,1)  ! INFO
1215      call interp_vert
1216     &    (told,var,lmold,llm,apsold,bpsold,aps,bps,
1217     &     psold,(imold+1)*(jmold+1))
1218      write (*,*) 'lect_start_archive: var ', var (1,jmold+1,1)  ! INFO
1219      call interp_horiz(var,t,imold,jmold,iim,jjm,llm,
1220     &                   rlonuold,rlatvold,rlonu,rlatv)
1221      write (*,*) 'lect_start_archive: t ', t(1,jjp1,1)  ! INFO
1222
1223c q2 : pbl wind variance
1224      write (*,*) 'lect_start_archive: q2old ', q2old (1,2,1)  ! INFO
1225      call interp_vert (q2old,varp1,lmold+1,llm+1,
1226     &     apsold,bpsold,aps,bps,psold,(imold+1)*(jmold+1))
1227      write (*,*) 'lect_start_archive: varp1 ', varp1 (1,2,1)  ! INFO
1228      call interp_horiz(varp1,q2s,imold,jmold,iim,jjm,llm+1,
1229     &                   rlonuold,rlatvold,rlonu,rlatv)
1230      write (*,*) 'lect_start_archive: q2s ', q2s (1,2,1)  ! INFO
1231      call gr_dyn_fi (llm+1,iim+1,jjm+1,ngridmx,q2s,q2)
1232      write (*,*) 'lect_start_archive: q2 ', q2 (1,2)  ! INFO
1233c     write(47,*) 'q2',q2
1234
1235c calcul des champ de vent; passage en vent covariant
1236      write (*,*) 'lect_start_archive: uold ', uold (1,2,1)  ! INFO
1237      call interp_vert
1238     & (uold,var,lmold,llm,apsold,bpsold,aps,bps,
1239     &  psold,(imold+1)*(jmold+1))
1240      write (*,*) 'lect_start_archive: var ', var (1,2,1)  ! INFO
1241      call interp_horiz(var,us,imold,jmold,iim,jjm,llm,
1242     &                   rlonuold,rlatvold,rlonu,rlatv)
1243      write (*,*) 'lect_start_archive: us ', us (1,2,1)   ! INFO
1244
1245      call interp_vert
1246     & (vold,var,lmold,llm,
1247     &  apsold,bpsold,aps,bps,psold,(imold+1)*(jmold+1))
1248      call interp_horiz(var,vs,imold,jmold,iim,jjm,llm,
1249     &                   rlonuold,rlatvold,rlonu,rlatv)
1250      call scal_wind(us,vs,unat,vnat)
1251      write (*,*) 'lect_start_archive: unat ', unat (1,2,1)    ! INFO
1252      do l=1,llm
1253        do j = 1, jjp1
1254          do i=1,iip1
1255            ucov( i,j,l ) = unat( i,j,l ) * cu(i,j)
1256c           ucov( i,j,l ) = 0
1257          end do
1258        end do
1259      end do
1260      write (*,*) 'lect_start_archive: ucov ', ucov (1,2,1)  ! INFO
1261c     write(48,*) 'ucov',ucov
1262      do l=1,llm
1263        do j = 1, jjm
1264          do i=1,iim
1265            vcov( i,j,l ) = vnat( i,j,l ) * cv(i,j)
1266c           vcov( i,j,l ) = 0
1267          end do
1268          vcov( iip1,j,l ) = vcov( 1,j,l )
1269        end do
1270      end do
1271c     write(49,*) 'ucov',vcov
1272
1273c traceurs surface
1274      do iq = 1, nqtot
1275            call interp_horiz(qsurfold(1,1,iq) ,qsurfs(1,1,iq),
1276     &                  imold,jmold,iim,jjm,1,
1277     &                  rlonuold,rlatvold,rlonu,rlatv)
1278      enddo
1279
1280      call gr_dyn_fi (nqtot,iim+1,jjm+1,ngridmx,qsurfs,qsurf)
1281
1282c traceurs 3D
1283      do  iq = 1, nqtot
1284            call interp_vert(qold(1,1,1,iq),var,lmold,llm,
1285     &        apsold,bpsold,aps,bps,psold,(imold+1)*(jmold+1))
1286            call interp_horiz(var,q(1,1,1,iq),imold,jmold,iim,jjm,llm,
1287     &                  rlonuold,rlatvold,rlonu,rlatv)
1288      enddo
1289cccccccccccccccccccccccccccccc     
1290c  make sure that sum of q = 1     
1291c dominent species is = 1 - sum(all other species)     
1292cccccccccccccccccccccccccccccc     
1293c      iqmax=1
1294c     
1295c      if (nqold.gt.10) then
1296c       do l=1,llm
1297c        do j=1,jjp1
1298c          do i=1,iip1
1299c           do iq=1,nqold
1300c            if (q(i,j,l,iq).gt.q(i,j,l,iqmax)) then
1301c              iqmax=iq
1302c            endif
1303c           enddo
1304c           q(i,j,l,iqmax)=1.
1305c           qtot(i,j,l)=0
1306c           do iq=1,nqold
1307c            if (iq.ne.iqmax) then       
1308c              q(i,j,l,iqmax)=q(i,j,l,iqmax)-q(i,j,l,iq)       
1309c            endif
1310c           enddo !iq
1311c           do iq=1,nqold
1312c            qtot(i,j,l)=qtot(i,j,l)+q(i,j,l,iq)
1313c            if (i.eq.1.and.j.eq.1.and.l.Eq.1) write(*,*)' qtot(i,j,l)',
1314c     $    qtot(i,j,l)
1315c           enddo !iq
1316c          enddo !i   
1317c         enddo !j   
1318c       enddo !l 
1319c      endif
1320ccccccccccccccccccccccccccccccc
1321
1322c     Periodicite :
1323      do  iq = 1, nqtot
1324         do l=1, llm
1325            do j = 1, jjp1
1326               q(iip1,j,l,iq) = q(1,j,l,iq)
1327            end do
1328         end do
1329      enddo
1330     
1331      call gr_dyn_fi (1,iim+1,jjm+1,ngridmx,co2ices,co2ice)
1332
1333c-----------------------------------------------------------------------
1334c   Initialisation  h:  (passage de t -> h)
1335c-----------------------------------------------------------------------
1336
1337      DO l=1,llm
1338         DO j=1,jjp1
1339            DO i=1,iim
1340               h(i,j,l) = t(i,j,l)*((ps(i,j)/preff)**kappa)
1341            END DO
1342            h(iip1,j,l) =  h(1,j,l)
1343         END DO
1344      END DO
1345
1346
1347c***********************************************************************
1348c***********************************************************************
1349c     Fin subroutine lecture ini
1350c***********************************************************************
1351c***********************************************************************
1352
1353      deallocate(timelist)
1354      deallocate(rlonuold, rlatvold)
1355      deallocate(rlonvold, rlatuold)
1356      deallocate(apsold,bpsold)
1357      deallocate(uold)
1358      deallocate(vold)
1359      deallocate(told)
1360      deallocate(psold)
1361      deallocate(phisold)
1362      deallocate(qold)
1363      deallocate(co2iceold)
1364      deallocate(tsurfold)
1365      deallocate(emisold)
1366      deallocate(q2old)
1367      deallocate(tsoilold)
1368      deallocate(tsoiloldnew)
1369      deallocate(inertiedatold)
1370      deallocate(inertiedatoldnew)
1371      deallocate(surfithold)
1372      deallocate(mlayerold)
1373      deallocate(qsurfold)
1374      deallocate(var,varp1)
1375
1376!      write(*,*)'lect_start_archive: END'
1377      return
1378      end
Note: See TracBrowser for help on using the repository browser.