source: trunk/LMDZ.PLUTO.old/libf/dyn3d/lect_start_archive.F @ 3436

Last change on this file since 3436 was 3175, checked in by emillour, 11 months ago

Pluto PCM:
Add the old Pluto LMDZ for reference (required prior step to making
an LMDZ.PLUTO using the same framework as the other physics packages).
TB+EM

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