Context Navigation

lect_start_archive.F @ 1670

Last change on this file since 1670 was 1647, checked in by jvatant, 8 years ago

+ Major clean of the new LMDZ.TITAN from too-generic options and routines (water, co2, ocean, surface type ...)
+ From this revision LMDZ.TITAN begins to be really separated from LMDZ.GENERIC
+ Partial desactivation of aerosols, only the dummy case is still enabled to keep the code running ( new aerosol routines to come in followings commits )

JVO

File size: 41.9 KB

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

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

Download in other formats:

Original Format