Context Navigation

lect_start_archive.F @ 1403

Last change on this file since 1403 was 1403, checked in by emillour, 10 years ago

All models: Reorganizing the physics/dynamics interface.

makelmdz and makelmdz_fcm scripts adapted to handle the new directory settings

misc: (replaces what was the "bibio" directory)
Should only contain extremely generic (and non physics or dynamics-specific) routines
Therefore moved initdynav.F90, initfluxsto.F, inithist.F, writedynav.F90, write_field.F90, writehist.F to "dyn3d_common"

dynlonlat_phylonlat: (new interface directory)
This directory contains routines relevent to physics/dynamics grid interactions, e.g. routines gr_dyn_fi or gr_fi_dyn and calfis
Moreover the dynlonlat_phylonlat contains directories "phy*" corresponding to each physics package "phy*" to be used. These subdirectories should only contain specific interfaces (e.g. iniphysiq) or main programs (e.g. newstart)

phy*/dyn1d: this subdirectory contains the 1D model using physics from phy*

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

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

Download in other formats:

Original Format