Context Navigation

lect_start_archive.F @ 222

Last change on this file since 222 was 174, checked in by emillour, 13 years ago

Mars GCM:

correct bug (introduced previously) in lect_start_archive.F on loop boundaries for soil temperature.

File size: 44.7 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	endif ! of if (depthinterpol)
526
527	c-----------------------------------------------------------------------
528	c 3.5 Read Soil thermal inertia
529	c-----------------------------------------------------------------------
530
531	ierr=NF_INQ_VARID(nid,"inertiedat",nvarid)
532	if (ierr.ne.NF_NOERR) then
533	write(,)'lect_start_archive: Could not find <inertiedat>'
534	write(,)' => Assuming this is an archive in old format'
535	therminertia_3D=.false.
536	write(,)' => Thermal inertia will be read from reference file'
537	volcapa=1.e6
538	write(,)' and soil volumetric heat capacity is set to ',
539	& volcapa
540	else
541	#ifdef NC_DOUBLE
542	ierr = NF_GET_VAR_DOUBLE(nid,nvarid,inertiedatold)
543	#else
544	ierr = NF_GET_VAR_REAL(nid,nvarid,inertiedatold)
545	#endif
546	if (ierr .NE. NF_NOERR) then
547	PRINT*, "lect_start_archive: Failed reading <inertiedat>"
548	CALL abort
549	endif
550	endif
551
552	c-----------------------------------------------------------------------
553	c 3.6 Lecture geopotentiel au sol
554	c-----------------------------------------------------------------------
555	c
556	ierr = NF_INQ_VARID (nid, "phisinit", nvarid)
557	IF (ierr .NE. NF_NOERR) THEN
558	PRINT*, "lect_start_archive: <phisinit> is missing"
559	CALL abort
560	ENDIF
561	#ifdef NC_DOUBLE
562	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, phisold)
563	#else
564	ierr = NF_GET_VAR_REAL(nid, nvarid, phisold)
565	#endif
566	IF (ierr .NE. NF_NOERR) THEN
567	PRINT*, "lect_start_archive: Failed loading <phisinit>"
568	CALL abort
569	ENDIF
570
571	C-----------------------------------------------------------------------
572	c lecture de "ptotalold" et "co2icetotalold"
573	c-----------------------------------------------------------------------
574	ptotalold = tab_cntrl(tab0+49)
575	co2icetotalold = tab_cntrl(tab0+50)
576
577	c-----------------------------------------------------------------------
578	c 4. Lecture du temps et choix
579	c-----------------------------------------------------------------------
580
581	c lecture du temps
582	c
583	ierr = NF_INQ_DIMID (nid, "Time", nvarid)
584	IF (ierr .NE. NF_NOERR) THEN
585	ierr = NF_INQ_DIMID (nid, "temps", nvarid)
586	IF (ierr .NE. NF_NOERR) THEN
587	PRINT*, "lect_start_archive: <Time> is missing"
588	CALL abort
589	endif
590	ENDIF
591
592	ierr = NF_INQ_VARID (nid, "Time", nvarid)
593	IF (ierr .NE. NF_NOERR) THEN
594	ierr = NF_INQ_VARID (nid, "temps", nvarid)
595	endif
596	#ifdef NC_DOUBLE
597	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, timelist)
598	#else
599	ierr = NF_GET_VAR_REAL(nid, nvarid, timelist)
600	#endif
601	IF (ierr .NE. NF_NOERR) THEN
602	PRINT*, "lect_start_archive: Failed loading <Time>"
603	CALL abort
604	ENDIF
605	c
606	write(,)
607	write(,)
608	write(,) 'Dates of the stored initial states:'
609	write(,) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
610	pi=2.*ASIN(1.)
611	do i=1,timelen
612	c call solarlong(timelist(i),sollong(i))
613	c sollong(i) = sollong(i)*180./pi
614	write(,) 'initial state at martian day: ',int(timelist(i))
615	c write(*,6) nint(timelist(i)),nint(mod(timelist(i),669)),
616	c . sollong(i)
617	end do
618
619	6 FORMAT(i7,i7,f9.3)
620
621	write(,)
622	write(,) 'Choose the martian day to use'
623	123 read(,,iostat=ierr) date
624	if(ierr.ne.0) goto 123
625	memo = 0
626	do i=1,timelen
627	if (date.eq.int(timelist(i))) then
628	memo = i
629	endif
630	end do
631
632	if (memo.eq.0) then
633	write(,)
634	write(,)
635	write(,) 'Wrong value for day number !!'
636	write(,)
637	write(,) 'Dates of the stored initial states:'
638	write(,) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
639	do i=1,timelen
640	write(,) 'initial state at martian day: ',nint(timelist(i))
641	c write(*,6) nint(timelist(i)),nint(mod(timelist(i),669))
642	end do
643	goto 123
644	endif
645
646	!-----------------------------------------------------------------------
647	! 5. Read (time-dependent) data from datafile
648	!-----------------------------------------------------------------------
649
650
651	c-----------------------------------------------------------------------
652	c 5.1 Lecture des champs 2D (co2ice, emis,ps,tsurf,Tg[10], q2surf)
653	c-----------------------------------------------------------------------
654
655	start=(/1,1,memo,0/)
656	count=(/imold+1,jmold+1,1,0/)
657
658	ierr = NF_INQ_VARID (nid, "co2ice", nvarid)
659	IF (ierr .NE. NF_NOERR) THEN
660	PRINT*, "lect_start_archive: <co2ice> is missing"
661	CALL abort
662	ENDIF
663	#ifdef NC_DOUBLE
664	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,co2iceold)
665	#else
666	ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,co2iceold)
667	#endif
668	IF (ierr .NE. NF_NOERR) THEN
669	PRINT*, "lect_start_archive: Failed loading <co2ice>"
670	PRINT*, NF_STRERROR(ierr)
671	CALL abort
672	ENDIF
673	c
674	ierr = NF_INQ_VARID (nid, "emis", nvarid)
675	IF (ierr .NE. NF_NOERR) THEN
676	PRINT*, "lect_start_archive: <emis> is missing"
677	CALL abort
678	ENDIF
679	#ifdef NC_DOUBLE
680	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,emisold)
681	#else
682	ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,emisold)
683	#endif
684	IF (ierr .NE. NF_NOERR) THEN
685	PRINT*, "lect_start_archive: Failed loading <emis>"
686	CALL abort
687	ENDIF
688	c
689	ierr = NF_INQ_VARID (nid, "ps", nvarid)
690	IF (ierr .NE. NF_NOERR) THEN
691	PRINT*, "lect_start_archive: <ps> is missing"
692	CALL abort
693	ENDIF
694	#ifdef NC_DOUBLE
695	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,psold)
696	#else
697	ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,psold)
698	#endif
699	IF (ierr .NE. NF_NOERR) THEN
700	PRINT*, "lect_start_archive: Failed loading <ps>"
701	CALL abort
702	ENDIF
703	c
704	ierr = NF_INQ_VARID (nid, "tsurf", nvarid)
705	IF (ierr .NE. NF_NOERR) THEN
706	PRINT*, "lect_start_archive: <tsurf> is missing"
707	CALL abort
708	ENDIF
709	#ifdef NC_DOUBLE
710	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,tsurfold)
711	#else
712	ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,tsurfold)
713	#endif
714	IF (ierr .NE. NF_NOERR) THEN
715	PRINT*, "lect_start_archive: Failed loading <tsurf>"
716	CALL abort
717	ENDIF
718	c
719	ierr = NF_INQ_VARID (nid, "q2surf", nvarid)
720	IF (ierr .NE. NF_NOERR) THEN
721	PRINT*, "lect_start_archive: <q2surf> is missing"
722	CALL abort
723	ENDIF
724	#ifdef NC_DOUBLE
725	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,q2old)
726	#else
727	ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,q2old)
728	#endif
729	IF (ierr .NE. NF_NOERR) THEN
730	PRINT*, "lect_start_archive: Failed loading <q2surf>"
731	CALL abort
732	ENDIF
733	c
734	write(,)"lect_start_archive: rlonuold:"
735	& ,rlonuold," rlatvold:",rlatvold
736	write(,)
737	c
738
739	c tracers: the 2 last ones are kept the 2 last one.
740	c the others keep their rank. ! No longer true.
741	c -------------------------------------------
742	! Surface tracers:
743	qsurfold(1:imold+1,1:jmold+1,1:nqmx)=0
744
745	DO iq=1,nqmx
746	IF (oldtracernames) THEN
747	txt=" "
748	write(txt,'(a5,i2.2)')'qsurf',iq
749	ELSE
750	txt=trim(tnom(iq))//"_surf"
751	if (txt.eq."h2o_vap") then
752	! There is no surface tracer for h2o_vap;
753	! "h2o_ice" should be loaded instead
754	txt="h2o_ice_surf"
755	write(,) 'lect_start_archive: loading surface tracer',
756	& ' h2o_ice instead of h2o_vap'
757	endif
758	ENDIF ! of IF (oldtracernames)
759	write(,) "lect_start_archive: loading tracer ",trim(txt)
760	ierr = NF_INQ_VARID (nid,txt,nvarid)
761	IF (ierr .NE. NF_NOERR) THEN
762	PRINT*, "lect_start_archive: ",
763	& " Tracer <",trim(txt),"> not found"
764	print*, "which (constant) value should it be initialized to?"
765	read(,) tmpval
766	qsurfold(1:imold+1,1:jmold+1,iq)=tmpval
767	ENDIF
768	#ifdef NC_DOUBLE
769	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,
770	& qsurfold(1,1,iq))
771	#else
772	ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,
773	& qsurfold(1,1,iq))
774	#endif
775	IF (ierr .NE. NF_NOERR) THEN
776	PRINT*, "lect_start_archive: ",
777	& " Failed loading <",trim(txt),">"
778	print*, "which (constant) value should it be initialized to?"
779	read(,) tmpval
780	qsurfold(1:imold+1,1:jmold+1,iq)=tmpval
781	ENDIF
782
783	ENDDO ! of DO iq=1,nqmx
784
785	!-----------------------------------------------------------------------
786	! 5.2 Read 3D subterranean fields
787	!-----------------------------------------------------------------------
788
789	start=(/1,1,1,memo/)
790	count=(/imold+1,jmold+1,nsoilold,1/)
791	!
792	! Read soil temperatures
793	!
794	if (olddepthdef) then ! tsoil stored using the 'old format'
795	start=(/1,1,memo,0/)
796	count=(/imold+1,jmold+1,1,0/) ! because the "Tg" are 2D datasets
797	do isoil=1,nsoilold
798	! write(,)'isoil',isoil
799	write(str2,'(i2.2)') isoil
800	c
801	ierr = NF_INQ_VARID (nid, "Tg"//str2, nvarid)
802	IF (ierr .NE. NF_NOERR) THEN
803	PRINT*, "lect_start_archive: ",
804	& "Field <","Tg"//str2,"> not found"
805	CALL abort
806	ENDIF
807	#ifdef NC_DOUBLE
808	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,
809	& tsoilold(1,1,isoil))
810	#else
811	ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,
812	& tsoilold(1,1,isoil))
813	#endif
814	IF (ierr .NE. NF_NOERR) THEN
815	PRINT*, "lect_start_archive: ",
816	& "Failed reading <","Tg"//str2,">"
817	CALL abort
818	ENDIF
819	c
820	enddo ! of do isoil=1,nsoilold
821
822	! reset 'start' and 'count' to "3D" behaviour
823	start=(/1,1,1,memo/)
824	count=(/imold+1,jmold+1,nsoilold,1/)
825
826	else
827	write(,) "lect_start_archive: loading tsoil "
828	ierr=NF_INQ_VARID(nid,"tsoil",nvarid)
829	if (ierr.ne.NF_NOERR) then
830	write(,)"lect_start_archive: Cannot find <tsoil>"
831	call abort
832	else
833	#ifdef NC_DOUBLE
834	ierr=NF_GET_VARA_DOUBLE(nid,nvarid,start,count,tsoilold)
835	#else
836	ierr=NF_GET_VARA_REAL(nid,nvarid,start,count,tsoilold)
837	#endif
838	endif ! of if (ierr.ne.NF_NOERR)
839
840	endif ! of if (olddepthdef)
841
842	!
843	! Read soil thermal inertias
844	!
845	! if (.not.olddepthdef) then ! no thermal inertia data in "old" archives
846	! ierr=NF_INQ_VARID(nid,"inertiedat",nvarid)
847	! if (ierr.ne.NF_NOERR) then
848	! write(,)"lect_start_archive: Cannot find <inertiedat>"
849	! call abort
850	! else
851	!#ifdef NC_DOUBLE
852	! ierr=NF_GET_VARA_DOUBLE(nid,nvarid,start,count,inertiedatold)
853	!#else
854	! ierr=NF_GET_VARA_REAL(nid,nvarid,start,count,inertiedatold)
855	!#endif
856	! endif ! of if (ierr.ne.NF_NOERR)
857	! endif
858
859	c-----------------------------------------------------------------------
860	c 5.3 Lecture des champs 3D (t,u,v, q2atm,q)
861	c-----------------------------------------------------------------------
862
863	start=(/1,1,1,memo/)
864	count=(/imold+1,jmold+1,lmold,1/)
865
866	c
867	ierr = NF_INQ_VARID (nid,"temp", nvarid)
868	IF (ierr .NE. NF_NOERR) THEN
869	PRINT*, "lect_start_archive: <temp> is missing"
870	CALL abort
871	ENDIF
872	#ifdef NC_DOUBLE
873	ierr = NF_GET_VARA_DOUBLE(nid, nvarid, start, count, told)
874	#else
875	ierr = NF_GET_VARA_REAL(nid, nvarid, start, count, told)
876	#endif
877	IF (ierr .NE. NF_NOERR) THEN
878	PRINT*, "lect_start_archive: Failed loading <temp>"
879	CALL abort
880	ENDIF
881	c
882	ierr = NF_INQ_VARID (nid,"u", nvarid)
883	IF (ierr .NE. NF_NOERR) THEN
884	PRINT*, "lect_start_archive: <u> is missing"
885	CALL abort
886	ENDIF
887	#ifdef NC_DOUBLE
888	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,uold)
889	#else
890	ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,uold)
891	#endif
892	IF (ierr .NE. NF_NOERR) THEN
893	PRINT*, "lect_start_archive: Failed loading <u>"
894	CALL abort
895	ENDIF
896	c
897	ierr = NF_INQ_VARID (nid,"v", nvarid)
898	IF (ierr .NE. NF_NOERR) THEN
899	PRINT*, "lect_start_archive: <v> is missing"
900	CALL abort
901	ENDIF
902	#ifdef NC_DOUBLE
903	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,vold)
904	#else
905	ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,vold)
906	#endif
907	IF (ierr .NE. NF_NOERR) THEN
908	PRINT*, "lect_start_archive: Failed loading <v>"
909	CALL abort
910	ENDIF
911	c
912	ierr = NF_INQ_VARID (nid,"q2atm", nvarid)
913	IF (ierr .NE. NF_NOERR) THEN
914	PRINT*, "lect_start_archive: <q2atm> is missing"
915	CALL abort
916	ENDIF
917	#ifdef NC_DOUBLE
918	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,q2old(1,1,2))
919	#else
920	ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,q2old(1,1,2))
921	#endif
922	IF (ierr .NE. NF_NOERR) THEN
923	PRINT*, "lect_start_archive: Failed loading <q2atm>"
924	CALL abort
925	ENDIF
926	c
927
928	c tracers: the 2 last ones are kept the 2 last one.
929	c the others keep their rank. ! No longer true.
930	c -------------------------------------------
931	! Tracers:
932	qold(1:imold+1,1:jmold+1,1:lmold,1:nqmx)=0
933
934	DO iq=1,nqmx
935	IF (oldtracernames) THEN
936	txt=" "
937	write(txt,'(a1,i2.2)')'q',iq
938	ELSE
939	txt=tnom(iq)
940	ENDIF
941	write(,)"lect_start_archive: loading tracer ",trim(txt)
942	ierr = NF_INQ_VARID (nid,txt,nvarid)
943	IF (ierr .NE. NF_NOERR) THEN
944	PRINT*, "lect_start_archive: ",
945	& " Tracer <",trim(txt),"> not found"
946	print*, "which (constant) value should it be initialized to?"
947	read(,) tmpval
948	qold(1:imold+1,1:jmold+1,1:lmold,iq)=tmpval
949	ENDIF
950	#ifdef NC_DOUBLE
951	ierr=NF_GET_VARA_DOUBLE(nid,nvarid,start,count,qold(1,1,1,iq))
952	#else
953	ierr=NF_GET_VARA_REAL(nid,nvarid,start,count,qold(1,1,1,iq))
954	#endif
955	IF (ierr .NE. NF_NOERR) THEN
956	PRINT*, "lect_start_archive: ",
957	& " Failed loading <",trim(txt),">"
958	print*, "which (constant) value should it be initialized to?"
959	read(,) tmpval
960	qold(1:imold+1,1:jmold+1,1:lmold,iq)=tmpval
961	ENDIF
962
963	ENDDO ! of DO iq=1,nqmx
964
965
966	c Chemin pour trouver les donnees de surface (albedo, relief, th.inertia...)
967	c -------------------------------------------------------------------------
968
969	! datapath = '/users/forget/gcm/data_mars_gcm'
970
971
972	!=======================================================================
973	! 6. Interpolation from old grid to new grid
974	!=======================================================================
975
976	c=======================================================================
977	c INTERPOLATION DANS LA NOUVELLE GRILLE et initialisation des variables
978	c=======================================================================
979	c Interpolation horizontale puis passage dans la grille physique pour
980	c les variables physique
981	c Interpolation verticale puis horizontale pour chaque variable 3D
982	c=======================================================================
983
984	c-----------------------------------------------------------------------
985	c 6.1 Variable 2d :
986	c-----------------------------------------------------------------------
987	c Relief
988	call interp_horiz (phisold,phisold_newgrid,imold,jmold,iim,jjm,1,
989	& rlonuold,rlatvold,rlonu,rlatv)
990
991	c Glace CO2
992	call interp_horiz (co2iceold,co2ices,imold,jmold,iim,jjm,1,
993	& rlonuold,rlatvold,rlonu,rlatv)
994
995	c Temperature de surface
996	call interp_horiz (tsurfold,tsurfs,imold,jmold,iim,jjm,1,
997	& rlonuold,rlatvold,rlonu,rlatv)
998	call gr_dyn_fi (1,iim+1,jjm+1,ngridmx,tsurfs,tsurf)
999	c write(44,*) 'tsurf', tsurf
1000
1001	c Temperature du sous-sol
1002	! call interp_horiz(tsoilold,tsoils,
1003	! & imold,jmold,iim,jjm,nsoilmx,
1004	! & rlonuold,rlatvold,rlonu,rlatv)
1005	! call gr_dyn_fi (nsoilmx,iim+1,jjm+1,ngridmx,tsoils,tsoil)
1006	c write(45,*) 'tsoil',tsoil
1007
1008	c Emissivite de la surface
1009	call interp_horiz (emisold,emiss,imold,jmold,iim,jjm,1,
1010	& rlonuold,rlatvold,rlonu,rlatv)
1011	call gr_dyn_fi (1,iim+1,jjm+1,ngridmx,emiss,emis)
1012	c write(46,*) 'emis',emis
1013	c-----------------------------------------------------------------------
1014	c 6.1.2 Traitement special de la pression au sol :
1015	c-----------------------------------------------------------------------
1016
1017	c Extrapolation la pression dans la nouvelle grille
1018	call interp_horiz(psold,ps,imold,jmold,iim,jjm,1,
1019	& rlonuold,rlatvold,rlonu,rlatv)
1020
1021	c-----------------------------------------------------------------------
1022	c On assure la conservation de la masse de l'atmosphere + calottes
1023	c-----------------------------------------------------------------------
1024
1025	ptotal = 0.
1026	co2icetotal = 0.
1027	DO j=1,jjp1
1028	DO i=1,iim
1029	ptotal=ptotal+ps(i,j)*aire(i,j)/g
1030	co2icetotal = co2icetotal + co2iceS(i,j)*aire(i,j)
1031	END DO
1032	END DO
1033
1034	write(,)
1035	write(,)'Old grid: mass of the atmosphere :',ptotalold
1036	write(,)'New grid: mass of the atmosphere :',ptotal
1037	write (,) 'Ratio new atm / old atm =', ptotal/ptotalold
1038	write(,)
1039	write(,)'Old grid: mass of CO2 ice:',co2icetotalold
1040	write(,)'New grid: mass of CO2 ice:',co2icetotal
1041	write(,)'Ratio new ice / old ice =',co2icetotal/co2icetotalold
1042	write(,)
1043
1044
1045	DO j=1,jjp1
1046	DO i=1,iip1
1047	ps(i,j)=ps(i,j) * ptotalold/ptotal
1048	END DO
1049	END DO
1050
1051	if ( co2icetotalold.gt.0.) then
1052	DO j=1,jjp1
1053	DO i=1,iip1
1054	co2iceS(i,j)=co2iceS(i,j) * co2icetotalold/co2icetotal
1055	END DO
1056	END DO
1057	end if
1058
1059	c-----------------------------------------------------------------------
1060	c 6.2 Subterranean 3d variables:
1061	c-----------------------------------------------------------------------
1062
1063	c-----------------------------------------------------------------------
1064	c 6.2.1 Thermal Inertia
1065	c Note: recall that inertiedat is a common in "comsoil.h"
1066	c-----------------------------------------------------------------------
1067
1068	! depth-wise interpolation, if required
1069	if (depthinterpol.and.(.not.olddepthdef)) then
1070	allocate(oldval(nsoilold))
1071	allocate(newval(nsoilmx))
1072	write(,)'lect_start_archive: WARNING: vertical interpolation o
1073	&f soil thermal inertia; might be wiser to reset it.'
1074	write(,)
1075
1076	do i=1,imold+1
1077	do j=1,jmold+1
1078	!copy old values
1079	oldval(1:nsoilold)=inertiedatold(i,j,1:nsoilold)
1080	!interpolate
1081	call interp_line(mlayerold,oldval,nsoilold,
1082	& mlayer,newval,nsoilmx)
1083	!copy interpolated values
1084	inertiedatoldnew(i,j,1:nsoilmx)=newval(1:nsoilmx)
1085	enddo
1086	enddo
1087	! cleanup
1088	deallocate(oldval)
1089	deallocate(newval)
1090	endif !of if (depthinterpol)
1091
1092	if (therminertia_3D) then
1093	! We have inertiedatold
1094	if((imold.ne.iim).or.(jmold.ne.jjm)) then
1095	write(,)'lect_start_archive: WARNING: horizontal interpolation
1096	&of thermal inertia; might be better to reset it.'
1097	write(,)
1098	endif
1099
1100	! Do horizontal interpolation
1101	if (depthinterpol) then
1102	call interp_horiz(inertiedatoldnew,inertiedatS,
1103	& imold,jmold,iim,jjm,nsoilmx,
1104	& rlonuold,rlatvold,rlonu,rlatv)
1105	else
1106	call interp_horiz(inertiedatold,inertiedatS,
1107	& imold,jmold,iim,jjm,nsoilold,
1108	& rlonuold,rlatvold,rlonu,rlatv)
1109	endif ! of if (depthinterpol)
1110
1111	else ! no 3D thermal inertia data
1112	write(,)'lect_start_archive: using reference surface inertia'
1113	! Use surface inertia (and extend it to all depths)
1114	do i=1,nsoilmx
1115	inertiedatS(1:iip1,1:jjp1,i)=surfith(1:iip1,1:jjp1)
1116	enddo
1117	! Build an old resolution surface thermal inertia
1118	! (will be needed for tsoil interpolation)
1119	call interp_horiz(surfith,surfithold,
1120	& iim,jjm,imold,jmold,1,
1121	& rlonu,rlatv,rlonuold,rlatvold)
1122	endif
1123
1124
1125	! Reshape inertiedatS to scalar grid as inertiedat
1126	call gr_dyn_fi (nsoilmx,iim+1,jjm+1,ngridmx,
1127	& inertiedatS,inertiedat)
1128
1129	c-----------------------------------------------------------------------
1130	c 6.2.2 Soil temperature
1131	c-----------------------------------------------------------------------
1132	! write(,) 'Soil'
1133	! Recast temperatures along soil depth, if necessary
1134	if (olddepthdef) then
1135	allocate(oldgrid(nsoilold+1))
1136	allocate(oldval(nsoilold+1))
1137	allocate(newval(nsoilmx))
1138	do i=1,imold+1
1139	do j=1,jmold+1
1140	! copy values
1141	oldval(1)=tsurfold(i,j)
1142	oldval(2:nsoilold+1)=tsoilold(i,j,1:nsoilold)
1143	! build vertical coordinate
1144	oldgrid(1)=0. ! ground
1145	oldgrid(2:nsoilold+1)=mlayerold(1:nsoilold)*
1146	& (surfithold(i,j)/1.e6)
1147	! Note; at this stage, we impose volcapa=1.e6 above
1148	! since volcapa isn't set in old soil definitions
1149
1150	! interpolate
1151	call interp_line(oldgrid,oldval,nsoilold+1,
1152	& mlayer,newval,nsoilmx)
1153	! copy result in tsoilold
1154	tsoiloldnew(i,j,1:nsoilmx)=newval(1:nsoilmx)
1155	enddo
1156	enddo
1157	! cleanup
1158	deallocate(oldgrid)
1159	deallocate(oldval)
1160	deallocate(newval)
1161
1162	else
1163	if (depthinterpol) then ! if vertical interpolation is required
1164	allocate(oldgrid(nsoilold+1))
1165	allocate(oldval(nsoilold+1))
1166	allocate(newval(nsoilmx))
1167	! build vertical coordinate
1168	oldgrid(1)=0. ! ground
1169	oldgrid(2:nsoilold+1)=mlayerold(1:nsoilold)
1170	do i=1,imold+1
1171	do j=1,jmold+1
1172	! copy values
1173	oldval(1)=tsurfold(i,j)
1174	oldval(2:nsoilold+1)=tsoilold(i,j,1:nsoilold)
1175	! interpolate
1176	call interp_line(oldgrid,oldval,nsoilold+1,
1177	& mlayer,newval,nsoilmx)
1178	! copy result in tsoilold
1179	tsoiloldnew(i,j,1:nsoilmx)=newval(1:nsoilmx)
1180	enddo
1181	enddo
1182	! write(,)'tsoiloldnew(1,1,1):',tsoiloldnew(1,1,1)
1183	! cleanup
1184	deallocate(oldgrid)
1185	deallocate(oldval)
1186	deallocate(newval)
1187
1188	else
1189	tsoiloldnew(:,:,:)=tsoilold(:,:,:)
1190	endif ! of if (depthinterpol)
1191	endif ! of if (olddepthdef)
1192
1193	! Do the horizontal interpolation
1194	call interp_horiz(tsoiloldnew,tsoilS,
1195	& imold,jmold,iim,jjm,nsoilmx,
1196	& rlonuold,rlatvold,rlonu,rlatv)
1197
1198	! Reshape tsoilS to scalar grid as tsoil
1199	call gr_dyn_fi (nsoilmx,iim+1,jjm+1,ngridmx,tsoilS,tsoil)
1200
1201
1202
1203	c-----------------------------------------------------------------------
1204	c 6.3 Variable 3d :
1205	c-----------------------------------------------------------------------
1206
1207	c temperatures atmospheriques
1208	write (,) 'lect_start_archive: told ', told (1,jmold+1,1) ! INFO
1209	call interp_vert
1210	& (told,var,lmold,llm,apsold,bpsold,aps,bps,
1211	& psold,(imold+1)*(jmold+1))
1212	write (,) 'lect_start_archive: var ', var (1,jmold+1,1) ! INFO
1213	call interp_horiz(var,t,imold,jmold,iim,jjm,llm,
1214	& rlonuold,rlatvold,rlonu,rlatv)
1215	write (,) 'lect_start_archive: t ', t(1,jjp1,1) ! INFO
1216
1217	c q2 : pbl wind variance
1218	write (,) 'lect_start_archive: q2old ', q2old (1,2,1) ! INFO
1219	call interp_vert (q2old,varp1,lmold+1,llm+1,
1220	& apsold,bpsold,aps,bps,psold,(imold+1)*(jmold+1))
1221	write (,) 'lect_start_archive: varp1 ', varp1 (1,2,1) ! INFO
1222	call interp_horiz(varp1,q2s,imold,jmold,iim,jjm,llm+1,
1223	& rlonuold,rlatvold,rlonu,rlatv)
1224	write (,) 'lect_start_archive: q2s ', q2s (1,2,1) ! INFO
1225	call gr_dyn_fi (llm+1,iim+1,jjm+1,ngridmx,q2s,q2)
1226	write (,) 'lect_start_archive: q2 ', q2 (1,2) ! INFO
1227	c write(47,*) 'q2',q2
1228
1229	c calcul des champ de vent; passage en vent covariant
1230	write (,) 'lect_start_archive: uold ', uold (1,2,1) ! INFO
1231	call interp_vert
1232	& (uold,var,lmold,llm,apsold,bpsold,aps,bps,
1233	& psold,(imold+1)*(jmold+1))
1234	write (,) 'lect_start_archive: var ', var (1,2,1) ! INFO
1235	call interp_horiz(var,us,imold,jmold,iim,jjm,llm,
1236	& rlonuold,rlatvold,rlonu,rlatv)
1237	write (,) 'lect_start_archive: us ', us (1,2,1) ! INFO
1238
1239	call interp_vert
1240	& (vold,var,lmold,llm,
1241	& apsold,bpsold,aps,bps,psold,(imold+1)*(jmold+1))
1242	call interp_horiz(var,vs,imold,jmold,iim,jjm,llm,
1243	& rlonuold,rlatvold,rlonu,rlatv)
1244	call scal_wind(us,vs,unat,vnat)
1245	write (,) 'lect_start_archive: unat ', unat (1,2,1) ! INFO
1246	do l=1,llm
1247	do j = 1, jjp1
1248	do i=1,iip1
1249	ucov( i,j,l ) = unat( i,j,l ) * cu(i,j)
1250	c ucov( i,j,l ) = 0
1251	end do
1252	end do
1253	end do
1254	write (,) 'lect_start_archive: ucov ', ucov (1,2,1) ! INFO
1255	c write(48,*) 'ucov',ucov
1256	do l=1,llm
1257	do j = 1, jjm
1258	do i=1,iim
1259	vcov( i,j,l ) = vnat( i,j,l ) * cv(i,j)
1260	c vcov( i,j,l ) = 0
1261	end do
1262	vcov( iip1,j,l ) = vcov( 1,j,l )
1263	end do
1264	end do
1265	c write(49,*) 'ucov',vcov
1266
1267	c traceurs surface
1268	do iq = 1, nqmx
1269	call interp_horiz(qsurfold(1,1,iq) ,qsurfs(1,1,iq),
1270	& imold,jmold,iim,jjm,1,
1271	& rlonuold,rlatvold,rlonu,rlatv)
1272	enddo
1273
1274	call gr_dyn_fi (nqmx,iim+1,jjm+1,ngridmx,qsurfs,qsurf)
1275
1276	c traceurs 3D
1277	do iq = 1, nqmx
1278	call interp_vert(qold(1,1,1,iq),var,lmold,llm,
1279	& apsold,bpsold,aps,bps,psold,(imold+1)*(jmold+1))
1280	call interp_horiz(var,q(1,1,1,iq),imold,jmold,iim,jjm,llm,
1281	& rlonuold,rlatvold,rlonu,rlatv)
1282	enddo
1283	cccccccccccccccccccccccccccccc
1284	c make sure that sum of q = 1
1285	c dominent species is = 1 - sum(all other species)
1286	cccccccccccccccccccccccccccccc
1287	c iqmax=1
1288	c
1289	c if (nqold.gt.10) then
1290	c do l=1,llm
1291	c do j=1,jjp1
1292	c do i=1,iip1
1293	c do iq=1,nqold
1294	c if (q(i,j,l,iq).gt.q(i,j,l,iqmax)) then
1295	c iqmax=iq
1296	c endif
1297	c enddo
1298	c q(i,j,l,iqmax)=1.
1299	c qtot(i,j,l)=0
1300	c do iq=1,nqold
1301	c if (iq.ne.iqmax) then
1302	c q(i,j,l,iqmax)=q(i,j,l,iqmax)-q(i,j,l,iq)
1303	c endif
1304	c enddo !iq
1305	c do iq=1,nqold
1306	c qtot(i,j,l)=qtot(i,j,l)+q(i,j,l,iq)
1307	c if (i.eq.1.and.j.eq.1.and.l.Eq.1) write(,)' qtot(i,j,l)',
1308	c $ qtot(i,j,l)
1309	c enddo !iq
1310	c enddo !i
1311	c enddo !j
1312	c enddo !l
1313	c endif
1314	ccccccccccccccccccccccccccccccc
1315
1316	c Periodicite :
1317	do iq = 1, nqmx
1318	do l=1, llm
1319	do j = 1, jjp1
1320	q(iip1,j,l,iq) = q(1,j,l,iq)
1321	end do
1322	end do
1323	enddo
1324
1325	call gr_dyn_fi (1,iim+1,jjm+1,ngridmx,co2ices,co2ice)
1326
1327	c-----------------------------------------------------------------------
1328	c Initialisation h: (passage de t -> h)
1329	c-----------------------------------------------------------------------
1330
1331	DO l=1,llm
1332	DO j=1,jjp1
1333	DO i=1,iim
1334	h(i,j,l) = t(i,j,l)((ps(i,j)/preff)*kappa)
1335	END DO
1336	h(iip1,j,l) = h(1,j,l)
1337	END DO
1338	END DO
1339
1340
1341	c***********************************************************************
1342	c***********************************************************************
1343	c Fin subroutine lecture ini
1344	c***********************************************************************
1345	c***********************************************************************
1346
1347	deallocate(timelist)
1348	deallocate(rlonuold, rlatvold)
1349	deallocate(rlonvold, rlatuold)
1350	deallocate(apsold,bpsold)
1351	deallocate(uold)
1352	deallocate(vold)
1353	deallocate(told)
1354	deallocate(psold)
1355	deallocate(phisold)
1356	deallocate(qold)
1357	deallocate(co2iceold)
1358	deallocate(tsurfold)
1359	deallocate(emisold)
1360	deallocate(q2old)
1361	deallocate(tsoilold)
1362	deallocate(tsoiloldnew)
1363	deallocate(inertiedatold)
1364	deallocate(inertiedatoldnew)
1365	deallocate(surfithold)
1366	deallocate(mlayerold)
1367	deallocate(qsurfold)
1368	deallocate(var,varp1)
1369
1370	! write(,)'lect_start_archive: END'
1371	return
1372	end

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

Download in other formats:

Original Format