Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

lect_start_archive.F @ 134

Last change on this file since 134 was 57, checked in by aslmd, 14 years ago
mineur LMD_MM_MARS: ajout du GCM ancienne physique, systeme maintenant complet sur SVN (ne manque que la base de donnees d'etats initiaux)
File size: 33.5 KB

Line
1
2	SUBROUTINE lect_start_archive(date,tsurf,tsoil,emis,q2,
3	. t,ucov,vcov,ps,co2ice,h,phisold_newgrid,q,qsurf,nid)
4	c=======================================================================
5	c
6	c
7	c Auteur: 05/1997 , 12/2003 : coord hybride FF
8	c ------
9	c
10	c
11	c Objet: Lecture des variables d'un fichier "start_archive"
12	c Plus besoin de régler ancienne valeurs grace
13	c a l'allocation dynamique de memoire (Yann Wanherdrick)
14	c
15	c
16	c
17	c=======================================================================
18
19	implicit none
20
21	#include "dimensions.h"
22	#include "dimphys.h"
23	#include "surfdat.h"
24	#include "dimradmars.h"
25	#include "yomaer.h"
26	#include "planete.h"
27	#include "paramet.h"
28	#include "comconst.h"
29	#include "comvert.h"
30	#include "comgeom2.h"
31	#include "control.h"
32	#include "logic.h"
33	#include "description.h"
34	#include "ener.h"
35	#include "temps.h"
36	#include "lmdstd.h"
37	#include "netcdf.inc"
38
39	c=======================================================================
40	c Declarations
41	c=======================================================================
42
43	c Variables dimension du fichier "ini"
44	c------------------------------------
45	INTEGER imold,jmold,lmold,nqold
46
47	c et autres:
48	c----------
49	INTEGER lnblnk
50	EXTERNAL lnblnk
51
52	c Variables pour les lectures des fichiers "ini"
53	c--------------------------------------------------
54	INTEGER sizei,timelen,dimid
55	INTEGER length
56	parameter (length = 100)
57	INTEGER tab0
58	INTEGER isoil,iq,iqmax
59	CHARACTER*2 str2
60
61	REAL dimfirst(4) ! tableau contenant les 1ers elements des dimensions
62
63	REAL dimlast(4) ! tableau contenant les derniers elements des dimensions
64
65	REAL dimcycl(4) ! tableau contenant les periodes des dimensions
66	CHARACTER*120 dimsource
67	CHARACTER*16 dimname
68	CHARACTER*80 dimtitle
69	CHARACTER*40 dimunits
70	INTEGER dimtype
71
72	INTEGER dimord(4) ! tableau contenant l''ordre
73	data dimord /1,2,3,4/ ! de sortie des dimensions
74
75	INTEGER vardim(4)
76	REAL date
77	INTEGER memo
78	character (len=50) :: tmpname
79
80	c Variable histoire
81	c------------------
82	REAL vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm) ! vents covariants
83	REAL h(iip1,jjp1,llm),ps(iip1,jjp1)
84	REAL q(iip1,jjp1,llm,nqmx),qtot(iip1,jjp1,llm)
85
86	c autre variables dynamique nouvelle grille
87	c------------------------------------------
88
89	c!-*-
90	integer klatdat,klongdat
91	PARAMETER (klatdat=180,klongdat=360)
92
93	c Physique sur grille scalaire
94	c----------------------------
95
96	c variable physique
97	c------------------
98	REAL tsurf(ngridmx),tsoil(ngridmx,nsoilmx),co2ice(ngridmx)
99	REAL emis(ngridmx)
100	REAL q2(ngridmx,nlayermx+1),qsurf(ngridmx,nqmx)
101	c REAL phisfi(ngridmx)
102
103	INTEGER i,j,l
104	INTEGER nid,nvarid
105	c REAL year_day,periheli,aphelie,peri_day
106	c REAL obliquit,z0,emin_turb,lmixmin
107	c REAL emissiv,emisice(2),albedice(2),tauvis
108	c REAL iceradius(2) , dtemisice(2)
109
110	EXTERNAL RAN1
111	REAL RAN1
112	EXTERNAL geopot,inigeom
113	integer ierr
114	integer ismin
115	external ismin
116	CHARACTER*80 datapath
117	integer, dimension(4) :: start,count
118
119	c Variable nouvelle grille naturelle au point scalaire
120	c------------------------------------------------------
121	real us(iip1,jjp1,llm),vs(iip1,jjp1,llm)
122	REAL phisold_newgrid(iip1,jjp1)
123	REAL t(iip1,jjp1,llm)
124	real tsurfS(iip1,jjp1),tsoilS(iip1,jjp1,nsoilmx)
125	real co2iceS(iip1,jjp1),emisS(iip1,jjp1)
126	REAL q2S(iip1,jjp1,llm+1),qsurfS(iip1,jjp1,nqmx)
127
128	real ptotal, co2icetotal
129
130	c Var intermediaires : vent naturel, mais pas coord scalaire
131	c-----------------------------------------------------------
132	real vnat(iip1,jjm,llm),unat(iip1,jjp1,llm)
133
134
135	c Variable de l'ancienne grille
136	c---------------------------------------------------------
137
138	real, dimension(:), allocatable :: timelist
139	real, dimension(:), allocatable :: rlonuold, rlatvold
140	real, dimension(:), allocatable :: rlonvold, rlatuold
141	real, dimension(:), allocatable :: apsold,bpsold
142	real, dimension(:,:,:), allocatable :: uold,vold,told,q2old
143	real, dimension(:,:,:), allocatable :: tsoilold,qsurfold
144	real, dimension(:,:), allocatable :: psold,phisold
145	real, dimension(:,:), allocatable :: co2iceold,tsurfold
146	real, dimension(:,:), allocatable :: emisold
147	real, dimension(:,:,:,:), allocatable :: qold
148
149	real tab_cntrl(100)
150
151	real ptotalold, co2icetotalold
152
153	c Variable intermediaires iutilise pour l'extrapolation verticale
154	c----------------------------------------------------------------
155	real, dimension(:,:,:), allocatable :: var,varp1
156
157	c=======================================================================
158
159	c Catching the axis lenghts for dynamic memory allocation
160
161	ierr= NF_INQ_DIMID(nid,"Time",dimid)
162	if (ierr.ne.NF_NOERR) then
163	ierr= NF_INQ_DIMID(nid,"temps",dimid)
164	endif
165	ierr= NF_INQ_DIMLEN(nid,dimid,timelen)
166
167	ierr= NF_INQ_DIMID(nid,"latitude",dimid)
168	if (ierr.ne.NF_NOERR) then
169	ierr= NF_INQ_DIMID(nid,"rlatu",dimid)
170	endif
171	ierr= NF_INQ_DIMLEN(nid,dimid,jmold)
172	jmold=jmold-1
173
174	ierr= NF_INQ_DIMID(nid,"longitude",dimid)
175	if (ierr.ne.NF_NOERR) then
176	ierr= NF_INQ_DIMID(nid,"rlonv",dimid)
177	endif
178	ierr= NF_INQ_DIMLEN(nid,dimid,imold)
179	imold=imold-1
180
181	ierr= NF_INQ_DIMID(nid,"altitude",dimid)
182	if (ierr.ne.NF_NOERR) then
183	ierr= NF_INQ_DIMID(nid,"sig_s",dimid)
184	endif
185	ierr= NF_INQ_DIMLEN(nid,dimid,lmold)
186
187	nqold=0
188	do
189	write(str2,'(i2.2)') nqold+1
190	ierr= NF_INQ_VARID(nid,'q'//str2,dimid)
191	! write(,) 'q'//str2
192	if (ierr.eq.NF_NOERR) then
193	nqold=nqold+1
194	else
195	exit
196	endif
197	enddo
198
199	write(,) "Start_archive dimensions:"
200	write(,) "longitude: ",imold
201	write(,) "latitude: ",jmold
202	write(,) "altitude: ",lmold
203	write(,) "tracers: ",nqold
204	write(,) "time lenght: ",timelen
205	write(,)
206
207	allocate(timelist(timelen))
208	allocate(rlonuold(imold+1), rlatvold(jmold))
209	allocate(rlonvold(imold+1), rlatuold(jmold+1))
210	allocate (apsold(lmold),bpsold(lmold))
211	allocate(uold(imold+1,jmold+1,lmold))
212	allocate(vold(imold+1,jmold+1,lmold))
213	allocate(told(imold+1,jmold+1,lmold))
214	allocate(psold(imold+1,jmold+1))
215	allocate(phisold(imold+1,jmold+1))
216	allocate(qold(imold+1,jmold+1,lmold,nqmx))
217	allocate(co2iceold(imold+1,jmold+1))
218	allocate(tsurfold(imold+1,jmold+1))
219	allocate(emisold(imold+1,jmold+1))
220	allocate(q2old(imold+1,jmold+1,lmold+1))
221	allocate(tsoilold(imold+1,jmold+1,nsoilmx))
222	allocate(qsurfold(imold+1,jmold+1,nqmx))
223
224	allocate(var (imold+1,jmold+1,llm))
225	allocate(varp1 (imold+1,jmold+1,llm+1))
226
227	write(,) 'q2',ngridmx,nlayermx+1
228	write(,) 'q2S',iip1,jjp1,llm+1
229	write(,) 'q2old',imold+1,jmold+1,lmold+1
230
231	C-----------------------------------------------------------------------
232	c Lecture du tableau des parametres du run
233	c (pour la lecture ulterieure de "ptotalold" et "co2icetotalold")
234	c-----------------------------------------------------------------------
235	c
236	ierr = NF_INQ_VARID (nid, "controle", nvarid)
237	IF (ierr .NE. NF_NOERR) THEN
238	PRINT*, "Lect_start_archive: champ <controle> est absent"
239	CALL abort
240	ENDIF
241	#ifdef NC_DOUBLE
242	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tab_cntrl)
243	#else
244	ierr = NF_GET_VAR_REAL(nid, nvarid, tab_cntrl)
245	#endif
246	IF (ierr .NE. NF_NOERR) THEN
247	PRINT*, "lect_start_archive: Lecture echoue pour <controle>"
248	CALL abort
249	ENDIF
250	c
251	tab0 = 50
252
253	c-----------------------------------------------------------------------
254	c Lecture des longitudes et latitudes
255	c-----------------------------------------------------------------------
256	c
257	ierr = NF_INQ_VARID (nid, "rlonv", nvarid)
258	IF (ierr .NE. NF_NOERR) THEN
259	PRINT*, "lect_start_archive: Le champ <rlonv> est absent"
260	CALL abort
261	ENDIF
262	#ifdef NC_DOUBLE
263	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlonvold)
264	#else
265	ierr = NF_GET_VAR_REAL(nid, nvarid, rlonvold)
266	#endif
267	IF (ierr .NE. NF_NOERR) THEN
268	PRINT*, "lect_start_archive: Lecture echouee pour <rlonv>"
269	CALL abort
270	ENDIF
271	c
272	ierr = NF_INQ_VARID (nid, "rlatu", nvarid)
273	IF (ierr .NE. NF_NOERR) THEN
274	PRINT*, "lect_start_archive: Le champ <rlatu> est absent"
275	CALL abort
276	ENDIF
277	#ifdef NC_DOUBLE
278	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlatuold)
279	#else
280	ierr = NF_GET_VAR_REAL(nid, nvarid, rlatuold)
281	#endif
282	IF (ierr .NE. NF_NOERR) THEN
283	PRINT*, "lect_start_archive: Lecture echouee pour <rlatu>"
284	CALL abort
285	ENDIF
286	c
287	ierr = NF_INQ_VARID (nid, "rlonu", nvarid)
288	IF (ierr .NE. NF_NOERR) THEN
289	PRINT*, "lect_start_archive: Le champ <rlonu> est absent"
290	CALL abort
291	ENDIF
292	#ifdef NC_DOUBLE
293	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlonuold)
294	#else
295	ierr = NF_GET_VAR_REAL(nid, nvarid, rlonuold)
296	#endif
297	IF (ierr .NE. NF_NOERR) THEN
298	PRINT*, "lect_start_archive: Lecture echouee pour <rlonu>"
299	CALL abort
300	ENDIF
301	c
302	ierr = NF_INQ_VARID (nid, "rlatv", nvarid)
303	IF (ierr .NE. NF_NOERR) THEN
304	PRINT*, "lect_start_archive: Le champ <rlatv> est absent"
305	CALL abort
306	ENDIF
307	#ifdef NC_DOUBLE
308	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlatvold)
309	#else
310	ierr = NF_GET_VAR_REAL(nid, nvarid, rlatvold)
311	#endif
312	IF (ierr .NE. NF_NOERR) THEN
313	PRINT*, "lect_start_archive: Lecture echouee pour <rlatv>"
314	CALL abort
315	ENDIF
316	c
317
318	c-----------------------------------------------------------------------
319	c Lecture des niveaux verticaux
320	c-----------------------------------------------------------------------
321	c
322	ierr = NF_INQ_VARID (nid, "aps", nvarid)
323	IF (ierr .NE. NF_NOERR) THEN
324	PRINT*, "lect_start_archive: Le champ <aps> est absent"
325	apsold=0
326	PRINT*, "<aps> set to 0"
327	ELSE
328	#ifdef NC_DOUBLE
329	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, apsold)
330	#else
331	ierr = NF_GET_VAR_REAL(nid, nvarid, apsold)
332	#endif
333	IF (ierr .NE. NF_NOERR) THEN
334	PRINT*, "lect_start_archive: Lecture echouee pour <aps>"
335	ENDIF
336	ENDIF
337	c
338	ierr = NF_INQ_VARID (nid, "bps", nvarid)
339	IF (ierr .NE. NF_NOERR) THEN
340	PRINT*, "lect_start_archive: Le champ <bps> est absent"
341	PRINT*, "It must be an old start_archive, lets look for sig_s"
342	ierr = NF_INQ_VARID (nid, "sig_s", nvarid)
343	IF (ierr .NE. NF_NOERR) THEN
344	PRINT*, "Nothing to do..."
345	CALL abort
346	ENDIF
347	ENDIF
348	#ifdef NC_DOUBLE
349	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, bpsold)
350	#else
351	ierr = NF_GET_VAR_REAL(nid, nvarid, bpsold)
352	#endif
353	IF (ierr .NE. NF_NOERR) THEN
354	PRINT*, "lect_start_archive: Lecture echouee pour <bps>"
355	CALL abort
356	END IF
357
358
359	c-----------------------------------------------------------------------
360	c Lecture geopotentiel au sol
361	c-----------------------------------------------------------------------
362	c
363	ierr = NF_INQ_VARID (nid, "phisinit", nvarid)
364	IF (ierr .NE. NF_NOERR) THEN
365	PRINT*, "lect_start_archive: Le champ <phisinit> est absent"
366	CALL abort
367	ENDIF
368	#ifdef NC_DOUBLE
369	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, phisold)
370	#else
371	ierr = NF_GET_VAR_REAL(nid, nvarid, phisold)
372	#endif
373	IF (ierr .NE. NF_NOERR) THEN
374	PRINT*, "lect_start_archive: Lecture echouee pour <phisinit>"
375	CALL abort
376	ENDIF
377
378	C-----------------------------------------------------------------------
379	c lecture de "ptotalold" et "co2icetotalold"
380	c-----------------------------------------------------------------------
381	ptotalold = tab_cntrl(tab0+49)
382	co2icetotalold = tab_cntrl(tab0+50)
383
384	c-----------------------------------------------------------------------
385	c Lecture du temps et choix
386	c-----------------------------------------------------------------------
387
388	c lecture du temps
389	c
390	ierr = NF_INQ_DIMID (nid, "Time", nvarid)
391	IF (ierr .NE. NF_NOERR) THEN
392	ierr = NF_INQ_DIMID (nid, "temps", nvarid)
393	IF (ierr .NE. NF_NOERR) THEN
394	PRINT*, "lect_start_archive: Le champ <Time> est absent"
395	CALL abort
396	endif
397	ENDIF
398
399	ierr = NF_INQ_VARID (nid, "Time", nvarid)
400	IF (ierr .NE. NF_NOERR) THEN
401	ierr = NF_INQ_VARID (nid, "temps", nvarid)
402	endif
403	#ifdef NC_DOUBLE
404	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, timelist)
405	#else
406	ierr = NF_GET_VAR_REAL(nid, nvarid, timelist)
407	#endif
408	IF (ierr .NE. NF_NOERR) THEN
409	PRINT*, "lect_start_archive: Lecture echouee pour <Time>"
410	CALL abort
411	ENDIF
412	c
413	write(,)
414	write(,)
415	write(,) 'Differentes dates des etats initiaux stockes:'
416	write(,) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
417	pi=2.*ASIN(1.)
418	do i=1,timelen
419	c call solarlong(timelist(i),sollong(i))
420	c sollong(i) = sollong(i)*180./pi
421	write(,) 'etat initial au jour martien' ,int(timelist(i))
422	c write(*,6) nint(timelist(i)),nint(mod(timelist(i),669)),
423	c . sollong(i)
424	end do
425
426	6 FORMAT(i7,i7,f9.3)
427
428	write(,)
429	write(,) 'Choix de la date'
430	123 read(,,iostat=ierr) date
431	if(ierr.ne.0) goto 123
432	memo = 0
433	do i=1,timelen
434	if (date.eq.int(timelist(i))) then
435	memo = i
436	endif
437	end do
438
439	if (memo.eq.0) then
440	write(,)
441	write(,)
442	write(,) 'He alors... Y sait pas lire !?!'
443	write(,)
444	write(,) 'Differentes dates des etats initiaux stockes:'
445	write(,) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
446	do i=1,timelen
447	write(,) 'etat initial au jour martien' ,nint(timelist(i))
448	c write(*,6) nint(timelist(i)),nint(mod(timelist(i),669))
449	end do
450	goto 123
451	endif
452
453	c-----------------------------------------------------------------------
454	c Lecture des champs 2D (co2ice, emis,ps,tsurf,Tg[10], q2surf)
455	c-----------------------------------------------------------------------
456
457	start=(/1,1,memo,0/)
458	count=(/imold+1,jmold+1,1,0/)
459
460	ierr = NF_INQ_VARID (nid, "co2ice", nvarid)
461	IF (ierr .NE. NF_NOERR) THEN
462	PRINT*, "lect_start_archive: Le champ <co2ice> est absent"
463	CALL abort
464	ENDIF
465	#ifdef NC_DOUBLE
466	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,co2iceold)
467	#else
468	ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,co2iceold)
469	#endif
470	IF (ierr .NE. NF_NOERR) THEN
471	PRINT*, "lect_start_archive: Lecture echouee pour <co2ice>"
472	PRINT*, NF_STRERROR(ierr)
473	CALL abort
474	ENDIF
475	c
476	ierr = NF_INQ_VARID (nid, "emis", nvarid)
477	IF (ierr .NE. NF_NOERR) THEN
478	PRINT*, "lect_start_archive: Le champ <emis> est absent"
479	CALL abort
480	ENDIF
481	#ifdef NC_DOUBLE
482	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,emisold)
483	#else
484	ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,emisold)
485	#endif
486	IF (ierr .NE. NF_NOERR) THEN
487	PRINT*, "lect_start_archive: Lecture echouee pour <emis>"
488	CALL abort
489	ENDIF
490	c
491	ierr = NF_INQ_VARID (nid, "ps", nvarid)
492	IF (ierr .NE. NF_NOERR) THEN
493	PRINT*, "lect_start_archive: Le champ <ps> est absent"
494	CALL abort
495	ENDIF
496	#ifdef NC_DOUBLE
497	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,psold)
498	#else
499	ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,psold)
500	#endif
501	IF (ierr .NE. NF_NOERR) THEN
502	PRINT*, "lect_start_archive: Lecture echouee pour <ps>"
503	CALL abort
504	ENDIF
505	c
506	ierr = NF_INQ_VARID (nid, "tsurf", nvarid)
507	IF (ierr .NE. NF_NOERR) THEN
508	PRINT*, "lect_start_archive: Le champ <tsurf> est absent"
509	CALL abort
510	ENDIF
511	#ifdef NC_DOUBLE
512	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,tsurfold)
513	#else
514	ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,tsurfold)
515	#endif
516	IF (ierr .NE. NF_NOERR) THEN
517	PRINT*, "lect_start_archive: Lecture echouee pour <tsurf>"
518	CALL abort
519	ENDIF
520	c
521	do isoil=1,nsoilmx
522	write(str2,'(i2.2)') isoil
523	c
524	ierr = NF_INQ_VARID (nid, "Tg"//str2, nvarid)
525	IF (ierr .NE. NF_NOERR) THEN
526	PRINT*, "lect_start_archive:
527	. Le champ <","Tg"//str2,"> est absent"
528	CALL abort
529	ENDIF
530	#ifdef NC_DOUBLE
531	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,
532	. tsoilold(1,1,isoil))
533	#else
534	ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,
535	. tsoilold(1,1,isoil))
536	#endif
537	IF (ierr .NE. NF_NOERR) THEN
538	PRINT*, "lect_start_archive:
539	. Lecture echouee pour <","Tg"//str2,">"
540	CALL abort
541	ENDIF
542	c
543	end do
544
545	c
546	ierr = NF_INQ_VARID (nid, "q2surf", nvarid)
547	IF (ierr .NE. NF_NOERR) THEN
548	PRINT*, "lect_start_archive: Le champ <q2surf> est absent"
549	CALL abort
550	ENDIF
551	#ifdef NC_DOUBLE
552	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,q2old)
553	#else
554	ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,q2old)
555	#endif
556	IF (ierr .NE. NF_NOERR) THEN
557	PRINT*, "lect_start_archive: Lecture echouee pour <q2surf>"
558	CALL abort
559	ENDIF
560	c
561	write (,) "rlonuold,rlatvold"
562	write (,) rlonuold
563	write (,) rlatvold
564	write (,)
565	c
566
567	c tracers: the 2 last ones are kept the 2 last one.
568	c the others keep their rank.
569	c -------------------------------------------
570
571	do iq=1,nqmx
572	call initial0((jmold+1)*(imold+1), qsurfold(1,1,iq))
573	enddo
574
575	iq=nqold
576	write(str2,'(i2.2)') iq
577	ierr = NF_INQ_VARID (nid, "qsurf"//str2, nvarid)
578	IF (ierr .NE. NF_NOERR) THEN
579	PRINT*, "lect_start_archive:
580	. Le champ <","qsurf"//str2,"> est absent"
581	CALL abort
582	ENDIF
583	#ifdef NC_DOUBLE
584	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,
585	. qsurfold(1,1,nqmx))
586	#else
587	ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,
588	. qsurfold(1,1,nqmx))
589	#endif
590	IF (ierr .NE. NF_NOERR) THEN
591	PRINT*, "lect_start_archive:
592	. Lecture echouee pour <","qsurf"//str2,">"
593	write (,) 'qsurf'//str2,' set to 0'
594	call initial0((jmold+1)*(imold+1), qsurfold(1,1,nqmx))
595	ENDIF
596
597	if ((nqold.gt.1).and.(nqmx.gt.1)) then
598	iq=nqold-1
599	write(str2,'(i2.2)') iq
600	ierr = NF_INQ_VARID (nid, "qsurf"//str2, nvarid)
601	IF (ierr .NE. NF_NOERR) THEN
602	PRINT*, "lect_start_archive:
603	. Le champ <","qsurf"//str2,"> est absent"
604	CALL abort
605	ENDIF
606	#ifdef NC_DOUBLE
607	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,
608	. qsurfold(1,1,nqmx-1))
609	#else
610	ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,
611	. qsurfold(1,1,nqmx-1))
612	#endif
613	IF (ierr .NE. NF_NOERR) THEN
614	PRINT*, "lect_start_archive:
615	. Lecture echouee pour <","qsurf"//str2,">"
616	write (,) 'qsurf'//str2,' set to 0'
617	call initial0((jmold+1)*(imold+1), qsurfold(1,1,nqmx-1))
618	ENDIF
619	endif
620
621	if (nqold.gt.2) then
622	do iq = 1, nqold-2
623	if (iq.lt.nqmx-1) then
624	write(str2,'(i2.2)') iq
625	ierr = NF_INQ_VARID (nid, "qsurf"//str2, nvarid)
626	IF (ierr .NE. NF_NOERR) THEN
627	PRINT*, "lect_start_archive:
628	. Le champ <","qsurf"//str2,"> est absent"
629	CALL abort
630	ENDIF
631	#ifdef NC_DOUBLE
632	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,
633	. qsurfold(1,1,iq))
634	#else
635	ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,
636	. qsurfold(1,1,iq))
637	#endif
638	IF (ierr .NE. NF_NOERR) THEN
639	PRINT*, "lect_start_archive:
640	. Lecture echouee pour <","qsurf"//str2,">"
641	write (,) 'qsurf'//str2,' set to 0'
642	call initial0((jmold+1)*(imold+1), qsurfold(1,1,iq))
643	ENDIF
644	end if
645	end do
646	end if
647
648	c-----------------------------------------------------------------------
649	c Lecture des champs 3D (t,u,v, q2atm,q)
650	c-----------------------------------------------------------------------
651
652	start=(/1,1,1,memo/)
653	count=(/imold+1,jmold+1,lmold,1/)
654
655	c
656	ierr = NF_INQ_VARID (nid,"temp", nvarid)
657	IF (ierr .NE. NF_NOERR) THEN
658	PRINT*, "lect_start_archive: Le champ <temp> est absent"
659	CALL abort
660	ENDIF
661	#ifdef NC_DOUBLE
662	ierr = NF_GET_VARA_DOUBLE(nid, nvarid, start, count, told)
663	#else
664	ierr = NF_GET_VARA_REAL(nid, nvarid, start, count, told)
665	#endif
666	IF (ierr .NE. NF_NOERR) THEN
667	PRINT*, "lect_start_archive: Lecture echouee pour <temp>"
668	CALL abort
669	ENDIF
670	c
671	ierr = NF_INQ_VARID (nid,"u", nvarid)
672	IF (ierr .NE. NF_NOERR) THEN
673	PRINT*, "lect_start_archive: Le champ <u> est absent"
674	CALL abort
675	ENDIF
676	#ifdef NC_DOUBLE
677	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,uold)
678	#else
679	ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,uold)
680	#endif
681	IF (ierr .NE. NF_NOERR) THEN
682	PRINT*, "lect_start_archive: Lecture echouee pour <u>"
683	CALL abort
684	ENDIF
685	c
686	ierr = NF_INQ_VARID (nid,"v", nvarid)
687	IF (ierr .NE. NF_NOERR) THEN
688	PRINT*, "lect_start_archive: Le champ <v> est absent"
689	CALL abort
690	ENDIF
691	#ifdef NC_DOUBLE
692	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,vold)
693	#else
694	ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,vold)
695	#endif
696	IF (ierr .NE. NF_NOERR) THEN
697	PRINT*, "lect_start_archive: Lecture echouee pour <v>"
698	CALL abort
699	ENDIF
700	c
701	ierr = NF_INQ_VARID (nid,"q2atm", nvarid)
702	IF (ierr .NE. NF_NOERR) THEN
703	PRINT*, "lect_start_archive: Le champ <q2atm> est absent"
704	CALL abort
705	ENDIF
706	#ifdef NC_DOUBLE
707	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,q2old(1,1,2))
708	#else
709	ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,q2old(1,1,2))
710	#endif
711	IF (ierr .NE. NF_NOERR) THEN
712	PRINT*, "lect_start_archive: Lecture echouee pour <q2atm>"
713	CALL abort
714	ENDIF
715	c
716
717	c tracers: the 2 last ones are kept the 2 last one.
718	c the others keep their rank.
719	c -------------------------------------------
720
721	do iq=1,nqmx
722	call initial0((jmold+1)(imold+1)lmold,qold(1,1,1,iq) )
723	enddo
724
725	iq=nqold
726	write(str2,'(i2.2)') iq
727	ierr = NF_INQ_VARID (nid, "q"//str2, nvarid)
728	IF (ierr .NE. NF_NOERR) THEN
729	PRINT*, "lect_start_archive:
730	. Le champ <","q"//str2,"> est absent"
731	CALL abort
732	ENDIF
733	#ifdef NC_DOUBLE
734	ierr= NF_GET_VARA_DOUBLE(nid,nvarid,start,count,qold(1,1,1,nqmx))
735	#else
736	ierr= NF_GET_VARA_REAL(nid,nvarid,start,count,qold(1,1,1,nqmx))
737	#endif
738	IF (ierr .NE. NF_NOERR) THEN
739	PRINT*, "lect_start_archive:
740	. Lecture echouee pour <","q"//str2,">"
741	write (,) 'q'//str2,' set to 1.E-30'
742	do l=1,lmold
743	do j=1,jmold+1
744	do i=1,imold+1
745	qold(1,1,1,nqmx)=1.e-30
746	end do
747	end do
748	end do
749
750	ENDIF
751
752	if ((nqold.gt.1).and.(nqmx.gt.1)) then
753	iq=nqold-1
754	write(str2,'(i2.2)') iq
755	ierr = NF_INQ_VARID (nid, "q"//str2, nvarid)
756	IF (ierr .NE. NF_NOERR) THEN
757	PRINT*, "lect_start_archive:
758	. Le champ <","q"//str2,"> est absent"
759	CALL abort
760	ENDIF
761	#ifdef NC_DOUBLE
762	ierr= NF_GET_VARA_DOUBLE(nid,nvarid,start,count,
763	. qold(1,1,1,nqmx-1))
764	#else
765	ierr= NF_GET_VARA_REAL(nid,nvarid,start,count,
766	. qold(1,1,1,nqmx-1))
767	#endif
768	IF (ierr .NE. NF_NOERR) THEN
769	PRINT*, "lect_start_archive:
770	. Lecture echouee pour <","q"//str2,">"
771	write (,) 'q'//str2,' set to 1.E-30'
772	do l=1,lmold
773	do j=1,jmold+1
774	do i=1,imold+1
775	qold(1,1,1,nqmx-1)=1.e-30
776	end do
777	end do
778	end do
779
780	ENDIF
781	endif
782
783	if (nqold.gt.2) then
784	do iq = 1, nqold-2
785	if (iq.lt.nqmx-1) then
786	write(str2,'(i2.2)') iq
787	ierr = NF_INQ_VARID (nid, "q"//str2, nvarid)
788	IF (ierr .NE. NF_NOERR) THEN
789	PRINT*, "lect_start_archive:
790	. Le champ <","q"//str2,"> est absent"
791	CALL abort
792	ENDIF
793	#ifdef NC_DOUBLE
794	ierr= NF_GET_VARA_DOUBLE(nid,nvarid,start,count,qold(1,1,1,iq))
795	#else
796	ierr= NF_GET_VARA_REAL(nid,nvarid,start,count,qold(1,1,1,iq))
797	#endif
798	IF (ierr .NE. NF_NOERR) THEN
799	PRINT*, "lect_start_archive:
800	. Lecture echouee pour <","q"//str2,">"
801	write (,) 'q'//str2,' set to 1.E-30 '
802	do l=1,lmold
803	do j=1,jmold+1
804	do i=1,imold+1
805	qold(1,1,1,iq)=1.e-30
806	end do
807	end do
808	end do
809
810	ENDIF
811	end if
812	end do
813	end if
814
815	c Chemin pour trouver les donnees de surface (albedo, relief, th.inertia...)
816	c -------------------------------------------------------------------------
817
818	datapath = '/users/forget/gcm/data_mars_gcm'
819
820
821	c=======================================================================
822	c INTERPOLATION DANS LA NOUVELLE GRILLE et initialisation des variables
823	c=======================================================================
824	c Interpolation horizontale puis passage dans la grille physique pour
825	c les variables physique
826	c Interpolation verticale puis horizontale pour chaque variable 3D
827	c=======================================================================
828
829	c-----------------------------------------------------------------------
830	c Variable 2d :
831	c-----------------------------------------------------------------------
832	c Relief
833	call interp_horiz (phisold,phisold_newgrid,imold,jmold,iim,jjm,1,
834	& rlonuold,rlatvold,rlonu,rlatv)
835
836	c Glace CO2
837	call interp_horiz (co2iceold,co2ices,imold,jmold,iim,jjm,1,
838	& rlonuold,rlatvold,rlonu,rlatv)
839
840	c Temperature de surface
841	call interp_horiz (tsurfold,tsurfs,imold,jmold,iim,jjm,1,
842	& rlonuold,rlatvold,rlonu,rlatv)
843	call gr_dyn_fi (1,iim+1,jjm+1,ngridmx,tsurfs,tsurf)
844	c write(44,*) 'tsurf', tsurf
845
846	c Temperature du sous-sol
847	call interp_horiz(tsoilold,tsoils,
848	& imold,jmold,iim,jjm,nsoilmx,
849	& rlonuold,rlatvold,rlonu,rlatv)
850	call gr_dyn_fi (nsoilmx,iim+1,jjm+1,ngridmx,tsoils,tsoil)
851	c write(45,*) 'tsoil',tsoil
852
853	c Emissivite de la surface
854	call interp_horiz (emisold,emiss,imold,jmold,iim,jjm,1,
855	& rlonuold,rlatvold,rlonu,rlatv)
856	call gr_dyn_fi (1,iim+1,jjm+1,ngridmx,emiss,emis)
857	c write(46,*) 'emis',emis
858	c-----------------------------------------------------------------------
859	c Traitement special de la pression au sol :
860	c-----------------------------------------------------------------------
861
862	c Extrapolation la pression dans la nouvelle grille
863	call interp_horiz(psold,ps,imold,jmold,iim,jjm,1,
864	& rlonuold,rlatvold,rlonu,rlatv)
865
866	c-----------------------------------------------------------------------
867	c On assure la conservation de la masse de l'atmosphere + calottes
868	c-----------------------------------------------------------------------
869
870	ptotal = 0.
871	co2icetotal = 0.
872	DO j=1,jjp1
873	DO i=1,iim
874	ptotal=ptotal+ps(i,j)*aire(i,j)/g
875	co2icetotal = co2icetotal + co2iceS(i,j)*aire(i,j)
876	END DO
877	END DO
878
879	write(,)
880	write(,)'Ancienne grille: masse de l atm :',ptotalold
881	write(,)'Nouvelle grille: masse de l atm :',ptotal
882	write (,) 'Ratio new atm./ old atm =', ptotal/ptotalold
883	write(,)
884	write(,)'Ancienne grille: masse de la glace CO2:',co2icetotalold
885	write(,)'Nouvelle grille: masse de la glace CO2:',co2icetotal
886	write(,)'Ratio new ice./old ice =',co2icetotal/co2icetotalold
887	write(,)
888
889
890	DO j=1,jjp1
891	DO i=1,iip1
892	ps(i,j)=ps(i,j) * ptotalold/ptotal
893	END DO
894	END DO
895
896	if ( co2icetotalold.gt.0.) then
897	DO j=1,jjp1
898	DO i=1,iip1
899	co2iceS(i,j)=co2iceS(i,j) * co2icetotalold/co2icetotal
900	END DO
901	END DO
902	end if
903
904	c-----------------------------------------------------------------------
905	c Variable 3d :
906	c-----------------------------------------------------------------------
907
908	c temperatures atmospheriques
909	write (,) 'told ', told (1,jmold+1,1) ! INFO
910	call interp_vert
911	& (told,var,lmold,llm,apsold,bpsold,aps,bps,
912	& psold,(imold+1)*(jmold+1))
913	write (,) 'var ', var (1,jmold+1,1) ! INFO
914	call interp_horiz(var,t,imold,jmold,iim,jjm,llm,
915	& rlonuold,rlatvold,rlonu,rlatv)
916	write (,) 't ', t(1,jjp1,1) ! INFO
917
918	c q2 : pbl wind variance
919	write (,) 'q2old ', q2old (1,2,1) ! INFO
920	call interp_vert (q2old,varp1,lmold+1,llm+1,
921	& apsold,bpsold,aps,bps,psold,(imold+1)*(jmold+1))
922	write (,) 'varp1 ', varp1 (1,2,1) ! INFO
923	call interp_horiz(varp1,q2s,imold,jmold,iim,jjm,llm+1,
924	& rlonuold,rlatvold,rlonu,rlatv)
925	write (,) 'q2s ', q2s (1,2,1) ! INFO
926	call gr_dyn_fi (llm+1,iim+1,jjm+1,ngridmx,q2s,q2)
927	write (,) 'q2 ', q2 (1,2) ! INFO
928	c write(47,*) 'q2',q2
929
930	c calcul des champ de vent; passage en vent covariant
931	write (,) 'uold ', uold (1,2,1) ! INFO
932	call interp_vert
933	& (uold,var,lmold,llm,apsold,bpsold,aps,bps,
934	& psold,(imold+1)*(jmold+1))
935	write (,) 'var ', var (1,2,1) ! INFO
936	call interp_horiz(var,us,imold,jmold,iim,jjm,llm,
937	& rlonuold,rlatvold,rlonu,rlatv)
938	write (,) 'us ', us (1,2,1) ! INFO
939
940	call interp_vert
941	& (vold,var,lmold,llm,
942	& apsold,bpsold,aps,bps,psold,(imold+1)*(jmold+1))
943	call interp_horiz(var,vs,imold,jmold,iim,jjm,llm,
944	& rlonuold,rlatvold,rlonu,rlatv)
945	call scal_wind(us,vs,unat,vnat)
946	write (,) 'unat ', unat (1,2,1) ! INFO
947	do l=1,llm
948	do j = 1, jjp1
949	do i=1,iip1
950	ucov( i,j,l ) = unat( i,j,l ) * cu(i,j)
951	c ucov( i,j,l ) = 0
952	end do
953	end do
954	end do
955	write (,) 'ucov ', ucov (1,2,1) ! INFO
956	c write(48,*) 'ucov',ucov
957	do l=1,llm
958	do j = 1, jjm
959	do i=1,iim
960	vcov( i,j,l ) = vnat( i,j,l ) * cv(i,j)
961	c vcov( i,j,l ) = 0
962	end do
963	vcov( iip1,j,l ) = vcov( 1,j,l )
964	end do
965	end do
966	c write(49,*) 'ucov',vcov
967
968	c traceurs surface
969	do iq = 1, nqmx
970	call interp_horiz(qsurfold(1,1,iq) ,qsurfs(1,1,iq),
971	& imold,jmold,iim,jjm,1,
972	& rlonuold,rlatvold,rlonu,rlatv)
973	enddo
974
975	call gr_dyn_fi (nqmx,iim+1,jjm+1,ngridmx,qsurfs,qsurf)
976
977	c traceurs 3D
978	do iq = 1, nqmx
979	call interp_vert(qold(1,1,1,iq),var,lmold,llm,
980	& apsold,bpsold,aps,bps,psold,(imold+1)*(jmold+1))
981	call interp_horiz(var,q(1,1,1,iq),imold,jmold,iim,jjm,llm,
982	& rlonuold,rlatvold,rlonu,rlatv)
983	enddo
984	cccccccccccccccccccccccccccccc
985	c make sure that sum of q = 1
986	c dominent species is = 1 - sum(all other species)
987	cccccccccccccccccccccccccccccc
988	c iqmax=1
989	c
990	c if (nqold.gt.10) then
991	c do l=1,llm
992	c do j=1,jjp1
993	c do i=1,iip1
994	c do iq=1,nqold
995	c if (q(i,j,l,iq).gt.q(i,j,l,iqmax)) then
996	c iqmax=iq
997	c endif
998	c enddo
999	c q(i,j,l,iqmax)=1.
1000	c qtot(i,j,l)=0
1001	c do iq=1,nqold
1002	c if (iq.ne.iqmax) then
1003	c q(i,j,l,iqmax)=q(i,j,l,iqmax)-q(i,j,l,iq)
1004	c endif
1005	c enddo !iq
1006	c do iq=1,nqold
1007	c qtot(i,j,l)=qtot(i,j,l)+q(i,j,l,iq)
1008	c if (i.eq.1.and.j.eq.1.and.l.Eq.1) write(,)' qtot(i,j,l)',
1009	c $ qtot(i,j,l)
1010	c enddo !iq
1011	c enddo !i
1012	c enddo !j
1013	c enddo !l
1014	c endif
1015	ccccccccccccccccccccccccccccccc
1016
1017	c Periodicite :
1018	do iq = 1, nqmx
1019	do l=1, llm
1020	do j = 1, jjp1
1021	q(iip1,j,l,iq) = q(1,j,l,iq)
1022	end do
1023	end do
1024	enddo
1025
1026	call gr_dyn_fi (1,iim+1,jjm+1,ngridmx,co2ices,co2ice)
1027
1028	c-----------------------------------------------------------------------
1029	c Initialisation h: (passage de t -> h)
1030	c-----------------------------------------------------------------------
1031
1032	DO l=1,llm
1033	DO j=1,jjp1
1034	DO i=1,iim
1035	h(i,j,l) = t(i,j,l)((ps(i,j)/preff)*kappa)
1036	END DO
1037	h(iip1,j,l) = h(1,j,l)
1038	END DO
1039	END DO
1040
1041
1042	c***********************************************************************
1043	c***********************************************************************
1044	c Fin subroutine lecture ini
1045	c***********************************************************************
1046	c***********************************************************************
1047
1048	deallocate(timelist)
1049	deallocate(rlonuold, rlatvold)
1050	deallocate(rlonvold, rlatuold)
1051	deallocate(apsold,bpsold)
1052	deallocate(uold)
1053	deallocate(vold)
1054	deallocate(told)
1055	deallocate(psold)
1056	deallocate(phisold)
1057	deallocate(qold)
1058	deallocate(co2iceold)
1059	deallocate(tsurfold)
1060	deallocate(emisold)
1061	deallocate(q2old)
1062	deallocate(tsoilold)
1063	deallocate(qsurfold)
1064	deallocate(var,varp1)
1065
1066	return
1067	end

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

Download in other formats:

Original Format