Context Navigation

lect_start_archive.F @ 2507

Last change on this file since 2507 was 2466, checked in by mvals, 4 years ago

Mars GCM:
Correction in lect_start_archive.F:

adds the variable "watercap" to insure its transmission in the case newstart.F uses start_archive.nc to create the start files.

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

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

Download in other formats:

Original Format