Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

lect_start_archive.F @ 91

Last change on this file since 91 was 38, checked in by emillour, 14 years ago
Ajout du modè Martien (mon LMDZ.MARS.BETA, du 28/01/2011) dans le rértoire mars, pour pouvoir suivre plus facilement les modifs. EM
File size: 48.9 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format