Context Navigation

lect_start_archive.F @ 536

Last change on this file since 536 was 253, checked in by emillour, 13 years ago

Generic GCM

Massive update to version 0.7

EM+RW

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

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

Download in other formats:

Original Format