Context Navigation

lect_start_archive.F @ 3026

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

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

Download in other formats:

Original Format