Context Navigation

lect_start_archive.F @ 1397

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

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

Download in other formats:

Original Format