Context Navigation

lect_start_archive.F @ 2322

Last change on this file since 2322 was 1711, checked in by mvals, 8 years ago

Mars GCM

implementation of a sub-grid water cloud fraction scheme (by A. Pottier): If CLFvarying is set to true in callphys.def, then the sub-grid cloud fraction routine is applied in watercloud.F and aeropacity.F.
accordingly modified files are: watercloud.F, aeropacity.F, callradite.F, conf_phys.F, phyetat0.F, phyredem.F90, physiq_mod.F, testphys1d.F, callkeys.h, newtsart.F, start2archive.F, lect_start_archive.F
added file: tcondwater.F90, used by watercloud.F to calculate the condensation temperature of water
watercloud.F, aeropacity.F, callradite.F are converted to module files as watercloud_mod.F, aeropacity_mod.F, callradite_mod.F

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

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

Download in other formats:

Original Format