Context Navigation

lect_start_archive.F @ 205

Last change on this file since 205 was 135, checked in by aslmd, 14 years ago
CHANGEMENT ARBORESCENCE ETAPE 2 -- NON COMPLET
File size: 43.8 KB

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

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

Download in other formats:

Original Format