Context Navigation

lect_start_archive.F @ 999

Last change on this file since 999 was 999, checked in by tnavarro, 11 years ago
Possibility to store multiple initial states in one start/startfi. This is RETROCOMPATIBLE. New option ecrithist in run.def to write data in start/startfi every ecrithist dynamical timestep. New option timestart in run.def to initialize the GCM with the time timestart stored in start
File size: 44.9 KB

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

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

Download in other formats:

Original Format