Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

lect_start_archive.F @ 3436

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

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

Download in other formats:

Original Format