Context Navigation

lect_start_archive.F @ 937

Last change on this file since 937 was 563, checked in by emillour, 13 years ago

Mars GCM:

Added option "q=profile" in newstart to initialize a given tracer with

a profile specified in an ASCII file.

Minor fix on lect_start_archive : when tracers were not found, a default

value was asked to user, but twice! Once is enough.

File size: 44.8 KB

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

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

Download in other formats:

Original Format