source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/lect_start_archive.F @ 224

Last change on this file since 224 was 222, checked in by ymipsl, 10 years ago

Creating temporary dynamico/lmdz/saturn branche

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