Context Navigation

lect_start_archive.F @ 3937

Last change on this file since 3937 was 3893, checked in by gmilcareck, 3 months ago
Remove all "call abort" and "stop" statement in LMDZ.GENERIC and replacing them by call abort_physic().
File size: 52.6 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	CALL abort_physic("lect_start_archive",
322	& "Field <controle> not found",1)
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	CALL abort_physic("lect_start_archive",
331	& "Lecture echoue pour <controle>",1)
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	CALL abort_physic("lect_start_archive",
343	& "Field <rlonv> not found",1)
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	CALL abort_physic("lect_start_archive",
352	& "Failed loading <rlonv>",1)
353	ENDIF
354	c
355	ierr = NF_INQ_VARID (nid, "rlatu", nvarid)
356	IF (ierr .NE. NF_NOERR) THEN
357	CALL abort_physic("lect_start_archive",
358	& "Field <rlatu> not found",1)
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	CALL abort_physic("lect_start_archive",
367	& "Failed loading <rlatu>",1)
368	ENDIF
369	c
370	ierr = NF_INQ_VARID (nid, "rlonu", nvarid)
371	IF (ierr .NE. NF_NOERR) THEN
372	CALL abort_physic("lect_start_archive",
373	& "Field <rlonu> not found",1)
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	CALL abort_physic("lect_start_archive",
382	& "Failed loading <rlonu>",1)
383	ENDIF
384	c
385	ierr = NF_INQ_VARID (nid, "rlatv", nvarid)
386	IF (ierr .NE. NF_NOERR) THEN
387	CALL abort_physic("lect_start_archive",
388	& "Field <rlatv> not found",1)
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	CALL abort_physic("lect_start_archive",
397	& "Failed loading <rlatv>",1)
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_physic("lect_start_archive",
429	& "Field <bps> or <sig_s> not found",1)
430	ENDIF
431	ENDIF
432	#ifdef NC_DOUBLE
433	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, bpsold)
434	#else
435	ierr = NF_GET_VAR_REAL(nid, nvarid, bpsold)
436	#endif
437	IF (ierr .NE. NF_NOERR) THEN
438	CALL abort_physic("lect_start_archive",
439	& "Failed loading <bps>",1)
440	END IF
441
442	c-----------------------------------------------------------------------
443	c 3.4 Read Soil layers depths
444	c-----------------------------------------------------------------------
445
446	ierr=NF_INQ_VARID(nid,"soildepth",nvarid)
447	if (ierr.ne.NF_NOERR) then
448	write(,)'lect_start_archive: Could not find <soildepth>'
449	write(,)' => Assuming this is an archive in old format'
450	olddepthdef=.true.
451	depthinterpol=.true.
452	! this is how soil depth was defined in ye old days
453	do isoil=1,nsoilold
454	mlayerold(isoil)=sqrt(887.75/3.14)((2.*(isoil-0.5))-1.)
455	enddo
456	else
457	#ifdef NC_DOUBLE
458	ierr = NF_GET_VAR_DOUBLE(nid,nvarid,mlayerold)
459	#else
460	ierr = NF_GET_VAR_REAL(nid,nvarid,mlayerold)
461	#endif
462	if (ierr .NE. NF_NOERR) then
463	CALL abort_physic("lect_start_archive",
464	& "Failed reading <soildepth>",1)
465	endif
466
467	endif !of if(ierr.ne.NF_NOERR)
468
469	! Read (or build) mlayer()
470	if (depthinterpol) then
471	! Build (default) new soil depths (mlayer(:) is in comsoil.h),
472	! as in soil_settings.F
473	write(,)' => Building default soil depths'
474	do isoil=0,nsoilmx-1
475	mlayer(isoil)=2.e-4(2.*(isoil-0.5))
476	enddo
477	write(,)' => mlayer: ',mlayer
478	! Also build (default) new soil interlayer depth layer(:)
479	do isoil=1,nsoilmx
480	layer(isoil)=sqrt(mlayer(0)mlayer(1))
481	& ((mlayer(1)/mlayer(0))**(isoil-1))
482	enddo
483	write(,)' => layer: ',layer
484	else ! read mlayer() from file
485	#ifdef NC_DOUBLE
486	ierr = NF_GET_VAR_DOUBLE(nid,nvarid,mlayer)
487	#else
488	ierr = NF_GET_VAR_REAL(nid,nvarid,mlayer)
489	#endif
490	if (ierr .NE. NF_NOERR) then
491	CALL abort_physic("lect_start_archive",
492	& "Failed reading <soildepth>",1)
493	endif
494	endif ! of if (depthinterpol)
495
496	c-----------------------------------------------------------------------
497	c 3.5 Read Soil thermal inertia
498	c-----------------------------------------------------------------------
499
500	ierr=NF_INQ_VARID(nid,"inertiedat",nvarid)
501	if (ierr.ne.NF_NOERR) then
502	write(,)'lect_start_archive: Could not find <inertiedat>'
503	write(,)' => Assuming this is an archive in old format'
504	therminertia_3D=.false.
505	write(,)' => Thermal inertia will be read from reference file'
506	volcapa=1.e6
507	write(,)' and soil volumetric heat capacity is set to ',
508	& volcapa
509	else
510	#ifdef NC_DOUBLE
511	ierr = NF_GET_VAR_DOUBLE(nid,nvarid,inertiedatold)
512	#else
513	ierr = NF_GET_VAR_REAL(nid,nvarid,inertiedatold)
514	#endif
515	if (ierr .NE. NF_NOERR) then
516	CALL abort_physic("lect_start_archive",
517	& "Failed reading <inertiedat>",1)
518	endif
519	endif
520
521	c-----------------------------------------------------------------------
522	c 3.6 Lecture geopotentiel au sol
523	c-----------------------------------------------------------------------
524	c
525	ierr = NF_INQ_VARID (nid, "phisinit", nvarid)
526	IF (ierr .NE. NF_NOERR) THEN
527	CALL abort_physic("lect_start_archive",
528	& "Field <phisinit> not found",1)
529	ENDIF
530	#ifdef NC_DOUBLE
531	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, phisold)
532	#else
533	ierr = NF_GET_VAR_REAL(nid, nvarid, phisold)
534	#endif
535	IF (ierr .NE. NF_NOERR) THEN
536	CALL abort_physic("lect_start_archive",
537	& "Failed loading <phisinit>",1)
538	ENDIF
539
540	C-----------------------------------------------------------------------
541	c lecture de "ptotalold" et "co2icetotalold"
542	c-----------------------------------------------------------------------
543	ptotalold = tab_cntrl(tab0+49)
544	co2icetotalold = tab_cntrl(tab0+50)
545
546	c-----------------------------------------------------------------------
547	c 4. Lecture du temps et choix
548	c-----------------------------------------------------------------------
549
550	c lecture du temps
551	c
552	ierr = NF_INQ_DIMID (nid, "Time", nvarid)
553	IF (ierr .NE. NF_NOERR) THEN
554	ierr = NF_INQ_DIMID (nid, "temps", nvarid)
555	IF (ierr .NE. NF_NOERR) THEN
556	CALL abort_physic("lect_start_archive",
557	& "Field <Time> not found",1)
558	endif
559	ENDIF
560
561	ierr = NF_INQ_VARID (nid, "Time", nvarid)
562	IF (ierr .NE. NF_NOERR) THEN
563	ierr = NF_INQ_VARID (nid, "temps", nvarid)
564	endif
565	#ifdef NC_DOUBLE
566	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, timelist)
567	#else
568	ierr = NF_GET_VAR_REAL(nid, nvarid, timelist)
569	#endif
570	IF (ierr .NE. NF_NOERR) THEN
571	CALL abort_physic("lect_start_archive",
572	& "Failed loading <Time>",1)
573	ENDIF
574	c
575	write(,)
576	write(,)
577	write(,) 'Available dates for the stored initial conditions:'
578	write(,) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
579	pi=2.*ASIN(1.)
580	do i=1,timelen
581	c call solarlong(timelist(i),sollong(i))
582	c sollong(i) = sollong(i)*180./pi
583	write(,) 'initial state for day ' ,int(timelist(i))
584	c write(*,6) nint(timelist(i)),nint(mod(timelist(i),669)),
585	c . sollong(i)
586	end do
587
588	6 FORMAT(i7,i7,f9.3)
589
590	write(,)
591	write(,) 'Choice for the date'
592	123 read(,,iostat=ierr) date
593	if(ierr.ne.0) goto 123
594	memo = 0
595	do i=1,timelen
596	if (date.eq.int(timelist(i))) then
597	memo = i
598	endif
599	end do
600
601	if (memo.eq.0) then
602	write(,)
603	write(,)
604	write(,) "Wrong value... can't you read !?!"
605	write(,)
606	write(,) 'Available dates for the stored initial conditions:'
607	write(,) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
608	do i=1,timelen
609	write(,) 'initial state for day ' ,nint(timelist(i))
610	c write(*,6) nint(timelist(i)),nint(mod(timelist(i),669))
611	end do
612	goto 123
613	endif
614
615	!-----------------------------------------------------------------------
616	! 5. Read (time-dependent) data from datafile
617	!-----------------------------------------------------------------------
618
619
620	c-----------------------------------------------------------------------
621	c 5.1 Lecture des champs 2D (co2ice, emis,ps,tsurf,Tg[10], qsurf)
622	c-----------------------------------------------------------------------
623
624	start=(/1,1,memo,0/)
625	count=(/imold+1,jmold+1,1,0/)
626
627	ierr = NF_INQ_VARID (nid, "emis", nvarid)
628	IF (ierr .NE. NF_NOERR) THEN
629	CALL abort_physic("lect_start_archive",
630	& "Field <emis> not found",1)
631	ENDIF
632	#ifdef NC_DOUBLE
633	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,emisold)
634	#else
635	ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,emisold)
636	#endif
637	IF (ierr .NE. NF_NOERR) THEN
638	CALL abort_physic("lect_start_archive",
639	& "Failed loading <emis>",1)
640	ENDIF
641	c
642	ierr = NF_INQ_VARID (nid, "ps", nvarid)
643	IF (ierr .NE. NF_NOERR) THEN
644	CALL abort_physic("lect_start_archive",
645	& "Field <ps> not found",1)
646	ENDIF
647	#ifdef NC_DOUBLE
648	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,psold)
649	#else
650	ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,psold)
651	#endif
652	IF (ierr .NE. NF_NOERR) THEN
653	CALL abort_physic("lect_start_archive",
654	& "Failed loading <ps>",1)
655	ENDIF
656	c
657	ierr = NF_INQ_VARID (nid, "tsurf", nvarid)
658	IF (ierr .NE. NF_NOERR) THEN
659	CALL abort_physic("lect_start_archive",
660	& "Field <tsurf> not found",1)
661	ENDIF
662	#ifdef NC_DOUBLE
663	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,tsurfold)
664	#else
665	ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,tsurfold)
666	#endif
667	IF (ierr .NE. NF_NOERR) THEN
668	CALL abort_physic("lect_start_archive",
669	& "Failed loading <tsurf>",1)
670	ENDIF
671	c
672	ierr = NF_INQ_VARID (nid, "q2surf", nvarid)
673	IF (ierr .NE. NF_NOERR) THEN
674	CALL abort_physic("lect_start_archive",
675	& "Field <q2surf> not found",1)
676	ENDIF
677	#ifdef NC_DOUBLE
678	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,q2old)
679	#else
680	ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,q2old)
681	#endif
682	IF (ierr .NE. NF_NOERR) THEN
683	CALL abort_physic("lect_start_archive",
684	& "Failed loading <q2surf>",1)
685	ENDIF
686	c
687	cc Slab ocean
688	if(ok_slab_ocean) then
689	start=(/1,1,1,memo/)
690	count=(/imold+1,jmold+1,nslay,1/)
691
692	ierr=NF_INQ_VARID(nid,"tslab",nvarid)
693	IF (ierr.ne.NF_NOERR) then
694	PRINT*,"lect_start_archive: Cannot find <tslab>"
695	ENDIF
696	#ifdef NC_DOUBLE
697	ierr=NF_GET_VARA_DOUBLE(nid,nvarid,start,count,tslabold)
698	#else
699	ierr=NF_GET_VARA_REAL(nid,nvarid,start,count,tslabold)
700	#endif
701	IF (ierr .NE. NF_NOERR) THEN
702	PRINT*, "lect_start_archive: Failed loading <tslab>"
703	ENDIF
704
705
706	c
707	start=(/1,1,memo,0/)
708	count=(/imold+1,jmold+1,1,0/)
709
710	ierr = NF_INQ_VARID (nid, "rnat", nvarid)
711	IF (ierr .NE. NF_NOERR) THEN
712	PRINT*, "lect_start_archive: Field <rnat> not found"
713	ENDIF
714	#ifdef NC_DOUBLE
715	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,rnatold)
716	#else
717	ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,rnatold)
718	#endif
719	IF (ierr .NE. NF_NOERR) THEN
720	PRINT*, "lect_start_archive: Failed loading <rnat>"
721	ENDIF
722	c
723	ierr = NF_INQ_VARID (nid, "pctsrf_sic", nvarid)
724	IF (ierr .NE. NF_NOERR) THEN
725	PRINT*, "lect_start_archive: Field <pctsrf_sic> not found"
726	ENDIF
727	#ifdef NC_DOUBLE
728	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,pctsrf_sicold)
729	#else
730	ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,pctsrf_sicold)
731	#endif
732	IF (ierr .NE. NF_NOERR) THEN
733	PRINT*, "lect_start_archive: Failed loading <pctsrf_sic>"
734	ENDIF
735	c
736	ierr = NF_INQ_VARID (nid, "tsea_ice", nvarid)
737	IF (ierr .NE. NF_NOERR) THEN
738	PRINT*, "lect_start_archive: Field <tsea_ice> not found"
739	ENDIF
740	#ifdef NC_DOUBLE
741	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,tsea_iceold)
742	#else
743	ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,tsea_iceold)
744	#endif
745	IF (ierr .NE. NF_NOERR) THEN
746	PRINT*, "lect_start_archive: Failed loading <tsea_ice>"
747	ENDIF
748	c
749	ierr = NF_INQ_VARID (nid, "tice", nvarid)
750	IF (ierr .NE. NF_NOERR) THEN
751	PRINT*, "lect_start_archive: Field <tice> not found"
752	ENDIF
753	#ifdef NC_DOUBLE
754	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,ticeold)
755	#else
756	ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,ticeold)
757	#endif
758	IF (ierr .NE. NF_NOERR) THEN
759	PRINT*, "lect_start_archive: Failed loading <tice>"
760	ENDIF
761	c
762	ierr = NF_INQ_VARID (nid, "sea_ice", nvarid)
763	IF (ierr .NE. NF_NOERR) THEN
764	PRINT*, "lect_start_archive: Field <sea_ice> not found"
765	ENDIF
766	#ifdef NC_DOUBLE
767	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,sea_iceold)
768	#else
769	ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,sea_iceold)
770	#endif
771	IF (ierr .NE. NF_NOERR) THEN
772	PRINT*, "lect_start_archive: Failed loading <sea_ice>"
773	ENDIF
774
775	ENDIF! ok_slab_ocean
776	c
777	write(,)"lect_start_archive: rlonuold:"
778	& ,rlonuold," rlatvold:",rlatvold
779	write(,)
780
781	! Surface tracers:
782	! initialize all surface tracers to zero
783	qsurfold(1:imold+1,1:jmold+1,1:nqtot)=0
784
785	DO iq=1,nqtot
786	txt=trim(tname(iq))//"_surf"
787	if (txt.eq."h2o_vap") then
788	! There is no surface tracer for h2o_vap;
789	! "h2o_ice" should be loaded instead
790	txt="h2o_ice_surf"
791	write(,) 'lect_start_archive: loading surface tracer',
792	& ' h2o_ice instead of h2o_vap'
793	endif
794
795
796	write(,) "lect_start_archive: loading tracer ",trim(txt)
797	ierr = NF_INQ_VARID (nid,txt,nvarid)
798	IF (ierr .NE. NF_NOERR) THEN
799	PRINT*, "lect_start_archive: ",
800	& " Tracer <",trim(txt),"> not found"
801
802	! print*,'RDW has added hack to let me continue...'
803	! CALL abort
804	ENDIF
805	#ifdef NC_DOUBLE
806	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,
807	& qsurfold(1,1,iq))
808	#else
809	ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,
810	& qsurfold(1,1,iq))
811	#endif
812	IF (ierr .NE. NF_NOERR) THEN
813	PRINT*, "lect_start_archive: ",
814	& " Failed loading <",trim(txt),">"
815	write (,) trim(txt),' is set to 0'
816	ENDIF
817
818	ENDDO ! of DO iq=1,nqtot
819
820
821	!-----------------------------------------------------------------------
822	! 5.2 Read 3D subterranean fields
823	!-----------------------------------------------------------------------
824
825	start=(/1,1,1,memo/)
826	count=(/imold+1,jmold+1,nsoilold,1/)
827	!
828	! Read soil temperatures
829	!
830	if (olddepthdef) then ! tsoil stored using the 'old format'
831	start=(/1,1,memo,0/)
832	count=(/imold+1,jmold+1,1,0/) ! because the "Tg" are 2D datasets
833	do isoil=1,nsoilold
834	write(str2,'(i2.2)') isoil
835	c
836	ierr = NF_INQ_VARID (nid, "Tg"//str2, nvarid)
837	IF (ierr .NE. NF_NOERR) THEN
838	PRINT*, "lect_start_archive: ",
839	& "Field <","Tg"//str2,"> not found"
840	CALL abort_physic("lect_start_archive",
841	& "Field <Tg> not found",1)
842	ENDIF
843	#ifdef NC_DOUBLE
844	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,
845	& tsoilold(1,1,isoil))
846	#else
847	ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,
848	& tsoilold(1,1,isoil))
849	#endif
850	IF (ierr .NE. NF_NOERR) THEN
851	PRINT*, "lect_start_archive: ",
852	& "Failed reading <","Tg"//str2,">"
853	CALL abort_physic("lect_start_archive",
854	& "Failed reading <Tg>",1)
855	ENDIF
856	c
857	enddo ! of do isoil=1,nsoilold
858
859	! reset 'start' and 'count' to "3D" behaviour
860	start=(/1,1,1,memo/)
861	count=(/imold+1,jmold+1,nsoilold,1/)
862
863	else
864	write(,) "lect_start_archive: loading tsoil "
865	ierr=NF_INQ_VARID(nid,"tsoil",nvarid)
866	if (ierr.ne.NF_NOERR) then
867	CALL abort_physic("lect_start_archive",
868	& "Cannot find <tsoil>",1)
869	else
870	#ifdef NC_DOUBLE
871	ierr=NF_GET_VARA_DOUBLE(nid,nvarid,start,count,tsoilold)
872	#else
873	ierr=NF_GET_VARA_REAL(nid,nvarid,start,count,tsoilold)
874	#endif
875	endif ! of if (ierr.ne.NF_NOERR)
876
877	endif ! of if (olddepthdef)
878
879	c-----------------------------------------------------------------------
880	c 5.3 Lecture des champs 3D (t,u,v, q2atm,q)
881	c-----------------------------------------------------------------------
882
883	start=(/1,1,1,memo/)
884	count=(/imold+1,jmold+1,lmold,1/)
885
886	c
887	ierr = NF_INQ_VARID (nid,"temp", nvarid)
888	IF (ierr .NE. NF_NOERR) THEN
889	CALL abort_physic("lect_start_archive",
890	& "Field <temp> not found",1)
891	ENDIF
892	#ifdef NC_DOUBLE
893	ierr = NF_GET_VARA_DOUBLE(nid, nvarid, start, count, told)
894	#else
895	ierr = NF_GET_VARA_REAL(nid, nvarid, start, count, told)
896	#endif
897	IF (ierr .NE. NF_NOERR) THEN
898	CALL abort_physic("lect_start_archive",
899	& "Failed loading <temp>",1)
900	ENDIF
901	c
902	ierr = NF_INQ_VARID (nid,"u", nvarid)
903	IF (ierr .NE. NF_NOERR) THEN
904	CALL abort_physic("lect_start_archive",
905	& "Field <u> not found",1)
906	ENDIF
907	#ifdef NC_DOUBLE
908	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,uold)
909	#else
910	ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,uold)
911	#endif
912	IF (ierr .NE. NF_NOERR) THEN
913	CALL abort_physic("lect_start_archive",
914	& "Failed loading <u>",1)
915	ENDIF
916	c
917	ierr = NF_INQ_VARID (nid,"v", nvarid)
918	IF (ierr .NE. NF_NOERR) THEN
919	CALL abort_physic("lect_start_archive",
920	& "Field <v> not found",1)
921	ENDIF
922	#ifdef NC_DOUBLE
923	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,vold)
924	#else
925	ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,vold)
926	#endif
927	IF (ierr .NE. NF_NOERR) THEN
928	CALL abort_physic("lect_start_archive",
929	& "Failed loading <v>",1)
930	ENDIF
931	c
932	ierr = NF_INQ_VARID (nid,"q2atm", nvarid)
933	IF (ierr .NE. NF_NOERR) THEN
934	CALL abort_physic("lect_start_archive",
935	& "Field <q2atm> not found",1)
936	ENDIF
937	#ifdef NC_DOUBLE
938	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,q2old(1,1,2))
939	#else
940	ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,q2old(1,1,2))
941	#endif
942	IF (ierr .NE. NF_NOERR) THEN
943	CALL abort_physic("lect_start_archive",
944	& "Failed loading <q2atm>",1)
945	ENDIF
946	c
947
948	! Tracers:
949	qold(1:imold+1,1:jmold+1,1:lmold,1:nqtot)=0.
950
951	DO iq=1,nqtot
952	txt=tname(iq)
953	write(,)"lect_start_archive: loading tracer ",trim(txt)
954	ierr = NF_INQ_VARID (nid,txt,nvarid)
955	IF (ierr .NE. NF_NOERR) THEN
956	PRINT*, "lect_start_archive: ",
957	& " Tracer <",trim(txt),"> not found"
958	! CALL abort
959	ENDIF
960	#ifdef NC_DOUBLE
961	ierr=NF_GET_VARA_DOUBLE(nid,nvarid,start,count,qold(1,1,1,iq))
962	#else
963	ierr=NF_GET_VARA_REAL(nid,nvarid,start,count,qold(1,1,1,iq))
964	#endif
965	IF (ierr .NE. NF_NOERR) THEN
966	PRINT*, "lect_start_archive: ",
967	& " Failed loading <",trim(txt),">"
968	write (,) trim(txt),' set to 1.E-30'
969	do l=1,lmold
970	do j=1,jmold+1
971	do i=1,imold+1
972	qold(i,j,l,iq)=1.e-30
973	end do
974	end do
975	end do
976	ENDIF
977
978	ENDDO ! of DO iq=1,nqtot
979
980	! Non-orographic GWs:
981	write(,)"lect_start_archive: loading du_nonoro_gwd"
982	ierr = NF_INQ_VARID (nid,"du_nonoro_gwd", nvarid)
983	IF (ierr .NE. NF_NOERR) THEN
984	PRINT*, "lect_start_archive: Field <du_nonoro_gwd> not found"
985	PRINT*, "Setting it to zero"
986	du_nonoro_gwdold(:,:,:)=0
987	ELSE
988	#ifdef NC_DOUBLE
989	ierr=NF_GET_VARA_DOUBLE(nid,nvarid,start,count,du_nonoro_gwdold)
990	#else
991	ierr=NF_GET_VARA_REAL(nid,nvarid,start,count,du_nonoro_gwdold)
992	#endif
993	IF (ierr .NE. NF_NOERR) THEN
994	CALL abort_physic("lect_start_archive",
995	& "Failed loading <du_nonoro_gwd>",1)
996	ENDIF
997	ENDIF
998
999	write(,)"lect_start_archive: loading dv_nonoro_gwd"
1000	ierr = NF_INQ_VARID (nid,"dv_nonoro_gwd", nvarid)
1001	IF (ierr .NE. NF_NOERR) THEN
1002	PRINT*, "lect_start_archive: Field <dv_nonoro_gwd> not found"
1003	PRINT*, "Setting it to zero"
1004	dv_nonoro_gwdold(:,:,:)=0
1005	ELSE
1006	#ifdef NC_DOUBLE
1007	ierr=NF_GET_VARA_DOUBLE(nid,nvarid,start,count,dv_nonoro_gwdold)
1008	#else
1009	ierr=NF_GET_VARA_REAL(nid,nvarid,start,count,dv_nonoro_gwdold)
1010	#endif
1011	IF (ierr .NE. NF_NOERR) THEN
1012	CALL abort_physic("lect_start_archive",
1013	& "Failed loading <dv_nonoro_gwd>",1)
1014	ENDIF
1015	ENDIF
1016
1017	write(,)"lect_start_archive: loading east_gwstress"
1018	ierr = NF_INQ_VARID (nid,"east_gwstress", nvarid)
1019	IF (ierr .NE. NF_NOERR) THEN
1020	PRINT*, "lect_start_archive: Field <east_gwstress> not found"
1021	PRINT*, "Setting it to zero"
1022	east_gwstressold(:,:,:)=0
1023	ELSE
1024	#ifdef NC_DOUBLE
1025	ierr=NF_GET_VARA_DOUBLE(nid,nvarid,start,count,east_gwstressold)
1026	#else
1027	ierr=NF_GET_VARA_REAL(nid,nvarid,start,count,east_gwstressold)
1028	#endif
1029	IF (ierr .NE. NF_NOERR) THEN
1030	CALL abort_physic("lect_start_archive",
1031	& "Failed loading <east_gwstress>",1)
1032	ENDIF
1033	ENDIF
1034
1035	write(,)"lect_start_archive: loading west_gwstress"
1036	ierr = NF_INQ_VARID (nid,"west_gwstress", nvarid)
1037	IF (ierr .NE. NF_NOERR) THEN
1038	PRINT*, "lect_start_archive: Field <west_gwstress> not found"
1039	PRINT*, "Setting it to zero"
1040	west_gwstressold(:,:,:)=0
1041	ELSE
1042	#ifdef NC_DOUBLE
1043	ierr=NF_GET_VARA_DOUBLE(nid,nvarid,start,count,west_gwstressold)
1044	#else
1045	ierr=NF_GET_VARA_REAL(nid,nvarid,start,count,west_gwstressold)
1046	#endif
1047	IF (ierr .NE. NF_NOERR) THEN
1048	CALL abort_physic("lect_start_archive",
1049	& "Failed loading <west_gwstress>",1)
1050	ENDIF
1051	ENDIF
1052
1053	!=======================================================================
1054	! 6. Interpolation from old grid to new grid
1055	!=======================================================================
1056
1057	c=======================================================================
1058	c INTERPOLATION DANS LA NOUVELLE GRILLE et initialisation des variables
1059	c=======================================================================
1060	c Interpolation horizontale puis passage dans la grille physique pour
1061	c les variables physique
1062	c Interpolation verticale puis horizontale pour chaque variable 3D
1063	c=======================================================================
1064
1065	c-----------------------------------------------------------------------
1066	c 6.1 Variable 2d :
1067	c-----------------------------------------------------------------------
1068	c Relief
1069	call interp_horiz (phisold,phisold_newgrid,imold,jmold,iim,jjm,1,
1070	& rlonuold,rlatvold,rlonu,rlatv)
1071
1072	! CO2 ice is now in qsurf(igcm_co2_ice)
1073	! call interp_horiz (co2iceold,co2ices,imold,jmold,iim,jjm,1,
1074	! & rlonuold,rlatvold,rlonu,rlatv)
1075
1076	c Temperature de surface
1077	call interp_horiz (tsurfold,tsurfs,imold,jmold,iim,jjm,1,
1078	& rlonuold,rlatvold,rlonu,rlatv)
1079	call gr_dyn_fi (1,iim+1,jjm+1,ngrid,tsurfs,tsurf)
1080	c write(44,*) 'tsurf', tsurf
1081
1082	c Temperature du sous-sol
1083	! call interp_horiz(tsoilold,tsoils,
1084	! & imold,jmold,iim,jjm,nsoilmx,
1085	! & rlonuold,rlatvold,rlonu,rlatv)
1086	! call gr_dyn_fi (nsoilmx,iim+1,jjm+1,ngrid,tsoils,tsoil)
1087	c write(45,*) 'tsoil',tsoil
1088
1089	c Emissivite de la surface
1090	call interp_horiz (emisold,emiss,imold,jmold,iim,jjm,1,
1091	& rlonuold,rlatvold,rlonu,rlatv)
1092	call gr_dyn_fi (1,iim+1,jjm+1,ngrid,emiss,emis)
1093	c write(46,*) 'emis',emis
1094
1095
1096
1097	c-----------------------------------------------------------------------
1098	c 6.1.2 Traitement special de la pression au sol :
1099	c-----------------------------------------------------------------------
1100
1101	c Extrapolation la pression dans la nouvelle grille
1102	call interp_horiz(psold,ps,imold,jmold,iim,jjm,1,
1103	& rlonuold,rlatvold,rlonu,rlatv)
1104
1105	c-----------------------------------------------------------------------
1106	c On assure la conservation de la masse de l'atmosphere + calottes
1107	c-----------------------------------------------------------------------
1108
1109	ptotal = 0.
1110	DO j=1,jjp1
1111	DO i=1,iim
1112	ptotal=ptotal+ps(i,j)*aire(i,j)/g
1113	END DO
1114	END DO
1115	co2icetotal = 0.
1116	if (igcm_co2_ice.ne.0) then
1117	! recast surface CO2 ice on new grid
1118	call interp_horiz(qsurfold(1,1,igcm_co2_ice),
1119	& qsurfs(1,1,igcm_co2_ice),
1120	& imold,jmold,iim,jjm,1,
1121	& rlonuold,rlatvold,rlonu,rlatv)
1122	DO j=1,jjp1
1123	DO i=1,iim
1124	!co2icetotal = co2icetotal + co2iceS(i,j)*aire(i,j)
1125	co2icetotal=co2icetotal+qsurfS(i,j,igcm_co2_ice)*aire(i,j)
1126	END DO
1127	END DO
1128	else
1129	write(,) "Warning: No co2_ice surface tracer"
1130	endif
1131
1132	write(,)
1133	write(,)'Old grid: atmospheric mass :',ptotalold
1134	write(,)'New grid: atmospheric mass :',ptotal
1135	write (,) 'Ratio new atm./ old atm =', ptotal/ptotalold
1136	write(,)
1137	write(,)'Old grid: mass of CO2 ice:',co2icetotalold
1138	write(,)'New grid: mass of CO2 ice:',co2icetotal
1139	if (co2icetotalold.ne.0.) then
1140	write(,)'Ratio new ice./old ice =',co2icetotal/co2icetotalold
1141	endif
1142	write(,)
1143
1144
1145	DO j=1,jjp1
1146	DO i=1,iip1
1147	ps(i,j)=ps(i,j) * ptotalold/ptotal
1148	END DO
1149	END DO
1150
1151	if ( co2icetotalold.gt.0.) then
1152	! DO j=1,jjp1
1153	! DO i=1,iip1
1154	! co2iceS(i,j)=co2iceS(i,j) * co2icetotalold/co2icetotal
1155	! END DO
1156	! END DO
1157	write(,) "Not enforcing conservation of surface CO2 ice"
1158	write(,) " (should be OK when CO2 ice is a tracer)"
1159	end if
1160
1161	c-----------------------------------------------------------------------
1162	c 6.2 Subterranean 3d variables:
1163	c-----------------------------------------------------------------------
1164
1165	c-----------------------------------------------------------------------
1166	c 6.2.1 Thermal Inertia
1167	c Note: recall that inertiedat is a common in "comsoil.h"
1168	c-----------------------------------------------------------------------
1169
1170	! depth-wise interpolation, if required
1171	if (depthinterpol.and.(.not.olddepthdef)) then
1172	allocate(oldval(nsoilold))
1173	allocate(newval(nsoilmx))
1174	write(,)'lect_start_archive: WARNING: vertical interpolation o
1175	&f soil thermal inertia; might be wiser to reset it.'
1176	write(,)
1177
1178	do i=1,imold+1
1179	do j=1,jmold+1
1180	!copy old values
1181	oldval(1:nsoilold)=inertiedatold(i,j,1:nsoilold)
1182	!interpolate
1183	call interp_line(mlayerold,oldval,nsoilold,
1184	& mlayer,newval,nsoilmx)
1185	!copy interpolated values
1186	inertiedatoldnew(i,j,1:nsoilmx)=newval(1:nsoilmx)
1187	enddo
1188	enddo
1189	! cleanup
1190	deallocate(oldval)
1191	deallocate(newval)
1192	endif !of if (depthinterpol)
1193
1194	if (therminertia_3D) then
1195	! We have inertiedatold
1196	if((imold.ne.iim).or.(jmold.ne.jjm)) then
1197	write(,)'lect_start_archive: WARNING: horizontal interpolation
1198	&of thermal inertia; might be better to reset it.'
1199	write(,)
1200	endif
1201
1202	! Do horizontal interpolation
1203	if (depthinterpol) then
1204	call interp_horiz(inertiedatoldnew,inertiedatS,
1205	& imold,jmold,iim,jjm,nsoilmx,
1206	& rlonuold,rlatvold,rlonu,rlatv)
1207	else
1208	call interp_horiz(inertiedatold,inertiedatS,
1209	& imold,jmold,iim,jjm,nsoilold,
1210	& rlonuold,rlatvold,rlonu,rlatv)
1211	endif ! of if (depthinterpol)
1212
1213	else ! no 3D thermal inertia data
1214	write(,)'lect_start_archive: using reference surface inertia'
1215	! Use surface inertia (and extend it to all depths)
1216	do i=1,nsoilmx
1217	inertiedatS(1:iip1,1:jjp1,i)=surfith(1:iip1,1:jjp1)
1218	enddo
1219	! Build an old resolution surface thermal inertia
1220	! (will be needed for tsoil interpolation)
1221	call interp_horiz(surfith,surfithold,
1222	& iim,jjm,imold,jmold,1,
1223	& rlonu,rlatv,rlonuold,rlatvold)
1224	endif
1225
1226
1227	! Reshape inertiedatS to scalar grid as inertiedat
1228	call gr_dyn_fi (nsoilmx,iim+1,jjm+1,ngrid,
1229	& inertiedatS,inertiedat)
1230
1231	c-----------------------------------------------------------------------
1232	c 6.2.2 Soil temperature
1233	c-----------------------------------------------------------------------
1234	! write(,) 'Soil'
1235
1236	!print*,'Problem in lect_start_archive interpolating'
1237	!print*,'to new resolution!!'
1238
1239	! Recast temperatures along soil depth, if necessary
1240	if (olddepthdef) then
1241	allocate(oldgrid(nsoilold+1))
1242	allocate(oldval(nsoilold+1))
1243	allocate(newval(nsoilmx))
1244	do i=1,imold+1
1245	do j=1,jmold+1
1246
1247	!if(i.gt.iip1 .or. j.gt.jjp1)then
1248	!print*,'Problem in lect_start_archive interpolating'
1249	!print*,'to new resolution!!'
1250	!call abort
1251	!endif
1252
1253	! copy values
1254	oldval(1)=tsurfold(i,j)
1255	! oldval(1)=tsurfS(i,j)
1256	oldval(2:nsoilold+1)=tsoilold(i,j,1:nsoilold)
1257	! build vertical coordinate
1258	oldgrid(1)=0. ! ground
1259	oldgrid(2:nsoilold+1)=mlayerold(1:nsoilold)*
1260	& (surfithold(i,j)/1.e6)
1261	! Note; at this stage, we impose volcapa=1.e6 above
1262	! since volcapa isn't set in old soil definitions
1263
1264	! interpolate
1265	call interp_line(oldgrid,oldval,nsoilold+1,
1266	& mlayer,newval,nsoilmx)
1267	! copy result in tsoilold
1268	tsoiloldnew(i,j,1:nsoilmx)=newval(1:nsoilmx)
1269	enddo
1270	enddo
1271	! cleanup
1272	deallocate(oldgrid)
1273	deallocate(oldval)
1274	deallocate(newval)
1275
1276	else
1277	if (depthinterpol) then ! if vertical interpolation is required
1278	allocate(oldgrid(nsoilold+1))
1279	allocate(oldval(nsoilold+1))
1280	allocate(newval(nsoilmx))
1281	! build vertical coordinate
1282	oldgrid(1)=0. ! ground
1283	oldgrid(2:nsoilold+1)=mlayerold(1:nsoilold)
1284	do i=1,imold+1
1285	do j=1,jmold+1
1286	! copy values
1287	oldval(1)=tsurfold(i,j)
1288	! oldval(1)=tsurfS(i,j)
1289	oldval(2:nsoilold+1)=tsoilold(i,j,1:nsoilold)
1290	! interpolate
1291	call interp_line(oldgrid,oldval,nsoilold+1,
1292	& mlayer,newval,nsoilmx)
1293	! copy result in tsoilold
1294	tsoiloldnew(i,j,1:nsoilmx)=newval(1:nsoilmx)
1295	enddo
1296	enddo
1297	! write(,)'tsoiloldnew(1,1,1):',tsoiloldnew(1,1,1)
1298	! cleanup
1299	deallocate(oldgrid)
1300	deallocate(oldval)
1301	deallocate(newval)
1302
1303	else
1304	tsoiloldnew(:,:,:)=tsoilold(:,:,:)
1305	endif ! of if (depthinterpol)
1306	endif ! of if (olddepthdef)
1307
1308	! Do the horizontal interpolation
1309	call interp_horiz(tsoiloldnew,tsoilS,
1310	& imold,jmold,iim,jjm,nsoilmx,
1311	& rlonuold,rlatvold,rlonu,rlatv)
1312
1313	! Reshape tsoilS to scalar grid as tsoil
1314	call gr_dyn_fi (nsoilmx,iim+1,jjm+1,ngrid,tsoilS,tsoil)
1315
1316	c-----------------------------------------------------------------------
1317	c 6.3 Slab Ocean :
1318	c-----------------------------------------------------------------------
1319	call interp_horiz (tslabold,tslabs,imold,jmold,iim,jjm,nslay,
1320	& rlonuold,rlatvold,rlonu,rlatv)
1321	call gr_dyn_fi (nslay,iim+1,jjm+1,ngrid,tslabs,tslab)
1322
1323	call interp_horiz (rnatold,rnats,imold,jmold,iim,jjm,1,
1324	& rlonuold,rlatvold,rlonu,rlatv)
1325	call gr_dyn_fi (1,iim+1,jjm+1,ngrid,rnats,rnat)
1326
1327	call interp_horiz (pctsrf_sicold,pctsrf_sics,imold,jmold,iim,
1328	& jjm,1,rlonuold,rlatvold,rlonu,rlatv)
1329	call gr_dyn_fi (1,iim+1,jjm+1,ngrid,pctsrf_sics,pctsrf_sic)
1330
1331	call interp_horiz (tsea_iceold,tsea_ices,imold,jmold,iim,jjm,1,
1332	& rlonuold,rlatvold,rlonu,rlatv)
1333	call gr_dyn_fi (1,iim+1,jjm+1,ngrid,tsea_ices,tsea_ice)
1334
1335	call interp_horiz (ticeold,tices,imold,jmold,iim,jjm,1,
1336	& rlonuold,rlatvold,rlonu,rlatv)
1337	call gr_dyn_fi (1,iim+1,jjm+1,ngrid,tices,tice)
1338
1339	call interp_horiz (sea_iceold,sea_ices,imold,jmold,iim,jjm,1,
1340	& rlonuold,rlatvold,rlonu,rlatv)
1341	call gr_dyn_fi (1,iim+1,jjm+1,ngrid,sea_ices,sea_ice)
1342
1343	c-----------------------------------------------------------------------
1344	c 6.4 Variable 3d :
1345	c-----------------------------------------------------------------------
1346
1347	c temperatures atmospheriques
1348	write (,) 'lect_start_archive: told ', told (1,jmold+1,1) ! INFO
1349	call interp_vert
1350	& (told,var,lmold,llm,apsold,bpsold,aps,bps,
1351	& psold,(imold+1)*(jmold+1))
1352	write (,) 'lect_start_archive: var ', var (1,jmold+1,1) ! INFO
1353	call interp_horiz(var,t,imold,jmold,iim,jjm,llm,
1354	& rlonuold,rlatvold,rlonu,rlatv)
1355	write (,) 'lect_start_archive: t ', t(1,jjp1,1) ! INFO
1356
1357	! Non-orographic GW
1358	call interp_vert
1359	& (du_nonoro_gwdold,var,lmold,llm,apsold,bpsold,aps,bps,
1360	& psold,(imold+1)*(jmold+1))
1361	call interp_horiz(var,du_nonoro_gwdS,imold,jmold,iim,jjm,llm,
1362	& rlonuold,rlatvold,rlonu,rlatv)
1363	call gr_dyn_fi(llm,iim+1,jjm+1,ngrid,du_nonoro_gwdS,du_nonoro_gwd)
1364
1365	call interp_vert
1366	& (dv_nonoro_gwdold,var,lmold,llm,apsold,bpsold,aps,bps,
1367	& psold,(imold+1)*(jmold+1))
1368	call interp_horiz(var,dv_nonoro_gwdS,imold,jmold,iim,jjm,llm,
1369	& rlonuold,rlatvold,rlonu,rlatv)
1370	call gr_dyn_fi(llm,iim+1,jjm+1,ngrid,dv_nonoro_gwdS,dv_nonoro_gwd)
1371
1372	call interp_vert
1373	& (east_gwstressold,var,lmold,llm,apsold,bpsold,aps,bps,
1374	& psold,(imold+1)*(jmold+1))
1375	call interp_horiz(var,east_gwstressS,imold,jmold,iim,jjm,llm,
1376	& rlonuold,rlatvold,rlonu,rlatv)
1377	call gr_dyn_fi(llm,iim+1,jjm+1,ngrid,east_gwstressS,east_gwstress)
1378
1379	call interp_vert
1380	& (west_gwstressold,var,lmold,llm,apsold,bpsold,aps,bps,
1381	& psold,(imold+1)*(jmold+1))
1382	call interp_horiz(var,west_gwstressS,imold,jmold,iim,jjm,llm,
1383	& rlonuold,rlatvold,rlonu,rlatv)
1384	call gr_dyn_fi(llm,iim+1,jjm+1,ngrid,west_gwstressS,west_gwstress)
1385
1386	c q2 : pbl wind variance
1387	write (,) 'lect_start_archive: q2old ', q2old (1,2,1) ! INFO
1388	call interp_vert (q2old,varp1,lmold+1,llm+1,
1389	& apsold,bpsold,ap,bp,psold,(imold+1)*(jmold+1))
1390	write (,) 'lect_start_archive: varp1 ', varp1 (1,2,1) ! INFO
1391	call interp_horiz(varp1,q2s,imold,jmold,iim,jjm,llm+1,
1392	& rlonuold,rlatvold,rlonu,rlatv)
1393	write (,) 'lect_start_archive: q2s ', q2s (1,2,1) ! INFO
1394	call gr_dyn_fi (llm+1,iim+1,jjm+1,ngrid,q2s,q2)
1395	write (,) 'lect_start_archive: q2 ', q2 (1,2) ! INFO
1396	c write(47,*) 'q2',q2
1397
1398	c calcul des champ de vent; passage en vent covariant
1399	write (,) 'lect_start_archive: uold ', uold (1,2,1) ! INFO
1400	call interp_vert
1401	& (uold,var,lmold,llm,apsold,bpsold,aps,bps,
1402	& psold,(imold+1)*(jmold+1))
1403	write (,) 'lect_start_archive: var ', var (1,2,1) ! INFO
1404	call interp_horiz(var,us,imold,jmold,iim,jjm,llm,
1405	& rlonuold,rlatvold,rlonu,rlatv)
1406	write (,) 'lect_start_archive: us ', us (1,2,1) ! INFO
1407
1408	call interp_vert
1409	& (vold,var,lmold,llm,
1410	& apsold,bpsold,aps,bps,psold,(imold+1)*(jmold+1))
1411	call interp_horiz(var,vs,imold,jmold,iim,jjm,llm,
1412	& rlonuold,rlatvold,rlonu,rlatv)
1413	call scal_wind(us,vs,unat,vnat)
1414	write (,) 'lect_start_archive: unat ', unat (1,2,1) ! INFO
1415	do l=1,llm
1416	do j = 1, jjp1
1417	do i=1,iip1
1418	ucov( i,j,l ) = unat( i,j,l ) * cu(i,j)
1419	c ucov( i,j,l ) = 0
1420	end do
1421	end do
1422	end do
1423	write (,) 'lect_start_archive: ucov ', ucov (1,2,1) ! INFO
1424	c write(48,*) 'ucov',ucov
1425	do l=1,llm
1426	do j = 1, jjm
1427	do i=1,iim
1428	vcov( i,j,l ) = vnat( i,j,l ) * cv(i,j)
1429	c vcov( i,j,l ) = 0
1430	end do
1431	vcov( iip1,j,l ) = vcov( 1,j,l )
1432	end do
1433	end do
1434	c write(49,*) 'ucov',vcov
1435
1436	if (nqtot .gt. 0) then
1437	c traceurs surface
1438	do iq = 1, nqtot
1439	call interp_horiz(qsurfold(1,1,iq) ,qsurfs(1,1,iq),
1440	& imold,jmold,iim,jjm,1,
1441	& rlonuold,rlatvold,rlonu,rlatv)
1442	enddo
1443
1444	call gr_dyn_fi (nqtot,iim+1,jjm+1,ngrid,qsurfs,qsurf)
1445
1446	c traceurs 3D
1447	do iq = 1, nqtot
1448	call interp_vert(qold(1,1,1,iq),var,lmold,llm,
1449	& apsold,bpsold,aps,bps,psold,(imold+1)*(jmold+1))
1450	call interp_horiz(var,q(1,1,1,iq),imold,jmold,iim,jjm,llm,
1451	& rlonuold,rlatvold,rlonu,rlatv)
1452	enddo
1453	cccccccccccccccccccccccccccccc
1454	c make sure that sum of q = 1
1455	c dominent species is = 1 - sum(all other species)
1456	cccccccccccccccccccccccccccccc
1457	c iqmax=1
1458	c
1459	c if (nqold.gt.10) then
1460	c do l=1,llm
1461	c do j=1,jjp1
1462	c do i=1,iip1
1463	c do iq=1,nqold
1464	c if (q(i,j,l,iq).gt.q(i,j,l,iqmax)) then
1465	c iqmax=iq
1466	c endif
1467	c enddo
1468	c q(i,j,l,iqmax)=1.
1469	c qtot(i,j,l)=0
1470	c do iq=1,nqold
1471	c if (iq.ne.iqmax) then
1472	c q(i,j,l,iqmax)=q(i,j,l,iqmax)-q(i,j,l,iq)
1473	c endif
1474	c enddo !iq
1475	c do iq=1,nqold
1476	c qtot(i,j,l)=qtot(i,j,l)+q(i,j,l,iq)
1477	c if (i.eq.1.and.j.eq.1.and.l.Eq.1) write(,)' qtot(i,j,l)',
1478	c $ qtot(i,j,l)
1479	c enddo !iq
1480	c enddo !i
1481	c enddo !j
1482	c enddo !l
1483	c endif
1484	ccccccccccccccccccccccccccccccc
1485
1486	c Periodicite :
1487	do iq = 1, nqtot
1488	do l=1, llm
1489	do j = 1, jjp1
1490	q(iip1,j,l,iq) = q(1,j,l,iq)
1491	end do
1492	end do
1493	enddo
1494
1495	! call gr_dyn_fi (1,iim+1,jjm+1,ngrid,co2ices,co2ice)
1496	! no need to transfer "co2ice" any more; it is in qsurf(igcm_co2_ice)
1497
1498	endif !! if nqtot .ne. 0
1499
1500	c-----------------------------------------------------------------------
1501	c Initialisation h: (passage de t -> h)
1502	c-----------------------------------------------------------------------
1503
1504	DO l=1,llm
1505	DO j=1,jjp1
1506	DO i=1,iim
1507	h(i,j,l) = t(i,j,l)((ps(i,j)/preff)*kappa)
1508	END DO
1509	h(iip1,j,l) = h(1,j,l)
1510	END DO
1511	END DO
1512
1513
1514	c***********************************************************************
1515	c***********************************************************************
1516	c Fin subroutine lecture ini
1517	c***********************************************************************
1518	c***********************************************************************
1519
1520	deallocate(timelist)
1521	deallocate(rlonuold, rlatvold)
1522	deallocate(rlonvold, rlatuold)
1523	deallocate(apsold,bpsold)
1524	deallocate(uold)
1525	deallocate(vold)
1526	deallocate(told)
1527	deallocate(psold)
1528	deallocate(phisold)
1529	deallocate(qold)
1530	deallocate(co2iceold)
1531	deallocate(tsurfold)
1532	deallocate(emisold)
1533	deallocate(q2old)
1534	deallocate(tsoilold)
1535	deallocate(tsoiloldnew)
1536	deallocate(inertiedatold)
1537	deallocate(inertiedatoldnew)
1538	deallocate(surfithold)
1539	deallocate(mlayerold)
1540	deallocate(qsurfold)
1541	deallocate(var,varp1)
1542	deallocate(tslabold)
1543	deallocate(rnatold)
1544	deallocate(pctsrf_sicold)
1545	deallocate(tsea_iceold)
1546	deallocate(ticeold)
1547	deallocate(sea_iceold)
1548
1549	end

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

Download in other formats:

Original Format