Context Navigation

lect_start_archive.F @ 1255

Last change on this file since 1255 was 1216, checked in by emillour, 11 years ago

Generic model:
Major cleanup, in order to ease the use of LMDZ.GENERIC with (parallel) dynamics
in LMDZ.COMMON: (NB: this will break LMDZ.UNIVERSAL, which should be thrashed
in the near future)

Updated makegcm_* scripts (and makdim) and added the "-full" (to enforce full recomputation of the model) option
In dyn3d: converted control.h to module control_mod.F90 and converted iniadvtrac.F to module infotrac.F90
Added module mod_const_mpi.F90 in dyn3d (not used in serial mode)
Rearanged input/outputs routines everywhere to handle serial/MPI cases. physdem.F => phyredem.F90 , phyetat0.F => phyetat0.F90 ; all read/write routines for startfi files are gathered in module iostart.F90
added parallelism related routines init_phys_lmdz.F90, comgeomphy.F90, dimphy.F90, iniphysiq.F90, mod_grid_phy_lmdz.F90, mod_phys_lmdz_mpi_data.F90, mod_phys_lmdz_mpi_transfert.F90, mod_phys_lmdz_omp_data.F90, mod_phys_lmdz_omp_transfert.F90, mod_phys_lmdz_para.F90, mod_phys_lmdz_transfert_para.F90 in phymars and mod_const_mpi.F90 in dyn3d (for compliance with parallelism)
added created generic routines 'planetwide_maxval' and 'planetwide_minval', in module "planetwide_mod", that enable obtaining the max and min of a field over the whole planet. This should be further imroved with computation of means (possibly area weighed), etc.

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

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

Download in other formats:

Original Format