Context Navigation

lect_start_archive.F @ 3452

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

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

Download in other formats:

Original Format