Context Navigation

lect_start_archive.F @ 1448

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

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

Download in other formats:

Original Format