Context Navigation

lect_start_archive.F @ 3555

Last change on this file since 3555 was 3423, checked in by bhatnags, 4 months ago
Generic-PCM Including variable "tice" in startarchive and newstart SB
File size: 52.1 KB

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

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

Download in other formats:

Original Format