Context Navigation

lect_start_archive.F @ 1130

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

Mars GCM:
Series of changes to enable running in parallel (using LMDZ.COMMON dynamics);
Current LMDZ.MARS can still notheless be compiled and run in serial mode
"as previously".
Summary of main changes:

Main programs (newstart, start2archive, xvik) that used to be in dyn3d have been moved to phymars.
dyn3d/control.h is now module control_mod.F90
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 parallel case)
created generic routines 'planetwide_maxval' and 'planetwide_minval', in module "planetwide_mod", that enable obtaining the min and max of a field over the whole planet.

File size: 45.3 KB

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

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

Download in other formats:

Original Format