source: trunk/LMDZ.GENERIC/libf/dyn3d/lect_start_archive.F @ 638

Last change on this file since 638 was 588, checked in by emillour, 13 years ago

Generic GCM:
Some cleanup and bug fixing:

  • "cloudfrac" was not well written to restartfi (wrong size).
  • missing save attribute for "reffrad" in physiq.F90.
  • cleanup recomputation of surface pressure in newstart and change loop order in interp_horiz (which "fixes" an odd behaviour which fills some arrays with zeros, but only when using some versions of ifort!)

EM

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