Context Navigation

lect_start_archive.F @ 1890

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

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

Download in other formats:

Original Format