Context Navigation

lect_start_archive.F @ 1237

Last change on this file since 1237 was 1232, checked in by emillour, 12 years ago
Mars GCM: Add "composition" option in newstart to enable global modification of CO2,N2,Ar,O2 and CO mixing ratios. Plus some cleanup on physics variables (saved arrays) which have been moved to modules. FF+EM
File size: 47.4 KB

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

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

Download in other formats:

Original Format