source: trunk/LMDZ.MARS/libf/phymars/lect_start_archive.F @ 1242

Last change on this file since 1242 was 1232, checked in by emillour, 11 years ago

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