Context Navigation

lect_start_archive.F @ 1303

Last change on this file since 1303 was 1297, checked in by bclmd, 10 years ago
LMDZ.GENERIC: slab ocean added
File size: 49.0 KB

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

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

Download in other formats:

Original Format