source: trunk/LMDZ.PLUTO/libf/dynphy_lonlat/phypluto/lect_start_archive.F @ 3546

Last change on this file since 3546 was 3275, checked in by afalco, 9 months ago

Pluto PCM:
Changed _vap to _gas;
Included surfprop.F90;
callcorrk includes methane
AF

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