Context Navigation

newstart.F @ 3556

Last change on this file since 3556 was 2405, checked in by slebonnois, 4 years ago
SL: correction of a bug in newstart, due to key cpofT missing ; + default is uniform T,u,v when vertical extension
File size: 42.3 KB

Line
1	C======================================================================
2	PROGRAM newstart
3	c=======================================================================
4	c
5	c
6	c Auteur: S. Lebonnois,
7	c a partir des newstart/start_archive/lect_start_archive martiens
8	c
9	c Derniere modif : 02/09 (ecriture des q*)
10	c 01/12 (inclusion dans svn dyn3d)
11	c
12	c Objet: Modify the grid for the initial state (LMD GCM VENUS/TITAN)
13	c ----- (from file NetCDF start_archive.nc)
14	c
15	c
16	c=======================================================================
17
18	use IOIPSL
19	USE filtreg_mod
20	USE startvar
21	USE control_mod
22	USE infotrac
23	use cpdet_mod, only: ini_cpdet,t2tpot
24	use exner_hyb_m, only: exner_hyb
25	use exner_milieu_m, only: exner_milieu
26	USE comconst_mod
27	USE comvert_mod, ONLY: ap,bp,presnivs,pa,preff,nivsigs,nivsig,
28	. aps,bps,scaleheight,pseudoalt,
29	. disvert_type,pressure_exner
30	USE serre_mod, ONLY: clon,clat,grossismx,grossismy,
31	& dzoomx,dzoomy,taux,tauy
32	USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
33	USE logic_mod, ONLY: iflag_trac,fxyhypb,ysinus
34	USE temps_mod, ONLY: day_ref,annee_ref
35	implicit none
36
37	#include "dimensions.h"
38	#include "paramet.h"
39	#include "comdissnew.h"
40	#include "comgeom2.h"
41	#include "description.h"
42	#include "dimsoil.h"
43	#include "netcdf.inc"
44
45	c-----------------------------------------------------------------------
46	c Declarations
47	c-----------------------------------------------------------------------
48
49	c Variables pour fichier "ini"
50	c------------------------------------
51	INTEGER imold,jmold,lmold,nqold,ip1jmp1old
52	INTEGER length
53	parameter (length = 100)
54	real tab_cntrl(2*length)
55	INTEGER isoil,iq,iqmax
56	CHARACTER*2 str2
57
58	c Variable histoire
59	c------------------
60	REAL vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm) ! vents covariants
61	REAL teta(iip1,jjp1,llm),ps(iip1,jjp1)
62	REAL phis(iip1,jjp1) ! geopotentiel au sol
63	REAL masse(ip1jmp1,llm) ! masse de l'atmosphere
64	REAL, ALLOCATABLE, DIMENSION(:,:,:,:):: q! champs advectes
65	REAL tab_cntrl_dyn(length) ! tableau des parametres de start
66
67	c variable physique
68	c------------------
69	integer ngridmx
70	parameter (ngridmx=(2+(jjm-1)*iim - 1/jjm))
71	REAL tab_cntrl_fi(length) ! tableau des parametres de startfi
72	real rlat(ngridmx),rlon(ngridmx)
73	REAL tsurf(ngridmx),tsoil(ngridmx,nsoilmx)
74	REAL albe(ngridmx),radsol(ngridmx),sollw(ngridmx)
75	real solsw(ngridmx),fder(ngridmx)
76	real sollwdown(ngridmx),dlw(ngridmx)
77	REAL zmea(ngridmx), zstd(ngridmx)
78	REAL zsig(ngridmx), zgam(ngridmx), zthe(ngridmx)
79	REAL zpic(ngridmx), zval(ngridmx)
80	real t_fi(ngridmx,llm)
81
82	c Variable nouvelle grille naturelle au point scalaire
83	c------------------------------------------------------
84	real us(iip1,jjp1,llm),vs(iip1,jjp1,llm)
85	REAL p3d(iip1,jjp1,llm+1) ! pression aux interfaces
86	REAL phisold_newgrid(iip1,jjp1)
87	REAL T(iip1,jjp1,llm)
88	real rlonS(iip1,jjp1),rlatS(iip1,jjp1)
89	real tsurfS(iip1,jjp1),tsoilS(iip1,jjp1,nsoilmx)
90	real albeS(ip1jmp1),radsolS(ip1jmp1),sollwS(ip1jmp1)
91	real solswS(ip1jmp1),fderS(ip1jmp1)
92	real sollwdownS(ip1jmp1),dlwS(ip1jmp1)
93	real zmeaS(ip1jmp1),zstdS(ip1jmp1),zsigS(ip1jmp1)
94	real zgamS(ip1jmp1),ztheS(ip1jmp1),zpicS(ip1jmp1)
95	real zvalS(ip1jmp1)
96
97	real ptotal
98
99	c Var intermediaires : vent naturel, mais pas coord scalaire
100	c-----------------------------------------------------------
101	real vnat(iip1,jjm,llm),unat(iip1,jjp1,llm)
102
103	REAL pks(iip1,jjp1) ! exner (f pour filtre)
104	REAL pk(iip1,jjp1,llm)
105	REAL pkf(iip1,jjp1,llm)
106	REAL alpha(iip1,jjp1,llm),beta(iip1,jjp1,llm)
107
108
109	c Variable de l'ancienne grille
110	c---------------------------------------------------------
111
112	real, dimension(:), allocatable :: rlonuold, rlatvold
113	real, dimension(:), allocatable :: rlonvold, rlatuold
114	real, dimension(:), allocatable :: nivsigsold,nivsigold
115	real, dimension(:), allocatable :: apold,bpold
116	real, dimension(:), allocatable :: presnivsold
117	real, dimension(:,:,:), allocatable :: uold,vold,told
118	real, dimension(:,:,:,:), allocatable :: qold
119	real, dimension(:,:,:), allocatable :: tsoilold
120	real, dimension(:,:), allocatable :: psold,phisold
121	real, dimension(:,:), allocatable :: tsurfold
122	real, dimension(:,:), allocatable :: albeold,radsolold
123	real, dimension(:,:), allocatable :: sollwold,solswold
124	real, dimension(:,:), allocatable :: fderold
125	real, dimension(:,:), allocatable :: dlwold,sollwdownold
126	real, dimension(:,:), allocatable :: zmeaold,zstdold,zsigold
127	real, dimension(:,:), allocatable :: zgamold,ztheold,zpicold
128	real, dimension(:,:), allocatable :: zvalold
129
130	real ptotalold
131
132	c Variable intermediaires iutilise pour l'extrapolation verticale
133	c----------------------------------------------------------------
134	real, dimension(:,:,:), allocatable :: var,varp1
135
136	c divers local
137	c-----------------
138
139	integer ierr,nid,nvarid
140	INTEGER ij, l,i,j
141	character*80 fichnom
142	integer, dimension(4) :: start,counter
143	REAL phisinverse(iip1,jjp1) ! geopotentiel au sol avant inversion
144	logical topoflag,notopo,albedoflag,razvitu,razvitv,uini
145	logical razTS,raztemp
146	real, dimension(:), allocatable :: tvira,dzst,zkm
147	real albedo
148
149	c=======================================================================
150	c INITIALISATIONS DIVERSES
151	c=======================================================================
152
153	c VENUS/TITAN
154
155	iflag_trac = 1
156	c-----------------------------------------------------------------------
157	c Initialisation des traceurs
158	c ---------------------------
159	c Choix du nombre de traceurs et du schema pour l'advection
160	c dans fichier traceur.def, par default ou via INCA
161	call infotrac_init
162
163	c Allocation de la tableau q : champs advectes
164	allocate(q(iip1,jjp1,llm,nqtot))
165
166	c-----------------------------------------------------------------------
167	c Ouverture du fichier a modifier (start_archive.nc)
168	c-----------------------------------------------------------------------
169
170	write(,) 'Creation d un etat initial a partir de'
171	write(,) './start_archive.nc'
172	write(,)
173	fichnom = 'start_archive.nc'
174	ierr = NF_OPEN (fichnom, NF_NOWRITE,nid)
175	IF (ierr.NE.NF_NOERR) THEN
176	write(6,*)' Pb d''ouverture du fichier ',fichnom
177	write(6,*)' ierr = ', ierr
178	CALL ABORT
179	ENDIF
180
181	c-----------------------------------------------------------------------
182	c Lecture du tableau des parametres du run (pour la dynamique)
183	c-----------------------------------------------------------------------
184
185	write(,) 'lecture tab_cntrl START_ARCHIVE'
186	c
187	ierr = NF_INQ_VARID (nid, "controle", nvarid)
188	#ifdef NC_DOUBLE
189	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tab_cntrl)
190	#else
191	ierr = NF_GET_VAR_REAL(nid, nvarid, tab_cntrl)
192	#endif
193	c
194	write(,) 'Impression de tab_cntrl'
195	do i=1,200
196	write(,) i,tab_cntrl(i)
197	enddo
198
199	c-----------------------------------------------------------------------
200	c Initialisation des constantes
201	c-----------------------------------------------------------------------
202
203	imold = tab_cntrl(1)
204	jmold = tab_cntrl(2)
205	lmold = tab_cntrl(3)
206	day_ref = tab_cntrl(4)
207	annee_ref = tab_cntrl(5)
208	rad = tab_cntrl(6)
209	omeg = tab_cntrl(7)
210	g = tab_cntrl(8)
211	cpp = tab_cntrl(9)
212	kappa = tab_cntrl(10)
213	daysec = tab_cntrl(11)
214	dtvr = tab_cntrl(12)
215	etot0 = tab_cntrl(13)
216	ptot0 = tab_cntrl(14)
217	ztot0 = tab_cntrl(15)
218	stot0 = tab_cntrl(16)
219	ang0 = tab_cntrl(17)
220	pa = tab_cntrl(18)
221	preff = tab_cntrl(19)
222	c
223	clon = tab_cntrl(20)
224	clat = tab_cntrl(21)
225	grossismx = tab_cntrl(22)
226	grossismy = tab_cntrl(23)
227
228	IF ( tab_cntrl(24).EQ.1. ) THEN
229	fxyhypb = . TRUE .
230	dzoomx = tab_cntrl(25)
231	dzoomy = tab_cntrl(26)
232	taux = tab_cntrl(28)
233	tauy = tab_cntrl(29)
234	ELSE
235	fxyhypb = . FALSE .
236	ysinus = . FALSE .
237	IF( tab_cntrl(27).EQ.1. ) ysinus = . TRUE.
238	ENDIF
239
240	ptotalold = tab_cntrl(2*length)
241
242	write(,) "Old dimensions:"
243	write(,) "longitude: ",imold
244	write(,) "latitude: ",jmold
245	write(,) "altitude: ",lmold
246	write(,)
247
248	ip1jmp1old = (imold+1-1/imold)*(jmold+1-1/jmold)
249
250	c dans run.def
251	CALL getin('planet_type',planet_type)
252	cpofT = .False.
253	if (planet_type.eq."venus") then
254	cpofT = .True.
255	endif
256	CALL getin('cpofT',cpofT)
257	call ini_cpdet
258
259	c=======================================================================
260	c CHANGEMENT DE CONSTANTES CONTENUES DANS tab_cntrl
261	c=======================================================================
262	c changement de la resolution dans le fichier de controle
263	tab_cntrl(1) = REAL(iim)
264	tab_cntrl(2) = REAL(jjm)
265	tab_cntrl(3) = REAL(llm)
266
267	c--------------------------------
268	c We use a specific run.def to read new parameters that need to be changed
269	c--------------------------------
270
271	c Changement de cpp:
272	CALL getin('cpp',cpp)
273	kappa = (8.314511/43.44E-3)/cpp
274	tab_cntrl(9) = cpp
275	tab_cntrl(10) = kappa
276
277	c CHANGING THE ZOOM PARAMETERS TO CHANGE THE GRID
278	CALL getin('clon',clon)
279	CALL getin('clat',clat)
280	tab_cntrl(20) = clon
281	tab_cntrl(21) = clat
282	CALL getin('grossismx',grossismx)
283	CALL getin('grossismy',grossismy)
284	tab_cntrl(22) = grossismx
285	tab_cntrl(23) = grossismy
286	CALL getin('fxyhypb',fxyhypb)
287	IF ( fxyhypb ) THEN
288	CALL getin('dzoomx',dzoomx)
289	CALL getin('dzoomy',dzoomy)
290	tab_cntrl(25) = dzoomx
291	tab_cntrl(26) = dzoomy
292	CALL getin('taux',taux)
293	CALL getin('tauy',tauy)
294	tab_cntrl(28) = taux
295	tab_cntrl(29) = tauy
296	ELSE
297	CALL getin('ysinus',ysinus)
298	IF (ysinus) THEN
299	tab_cntrl(27) = 1
300	ELSE
301	tab_cntrl(27) = 0
302	ENDIF
303	ENDIF
304
305	c changes must be done BEFORE these lines...
306	DO l=1,length
307	tab_cntrl_dyn(l)= tab_cntrl(l)
308	tab_cntrl_fi(l) = tab_cntrl(l+length)
309	ENDDO
310
311	c-----------------------------------------------------------------------
312	c Autres initialisations
313	c-----------------------------------------------------------------------
314
315	CALL iniconst
316	CALL inigeom
317	call inifilr
318
319	c-----------------------------------------------------------------------
320	c Allocations des anciennes variables
321	c-----------------------------------------------------------------------
322
323	allocate(rlonuold(imold+1), rlatvold(jmold))
324	allocate(rlonvold(imold+1), rlatuold(jmold+1))
325	allocate(nivsigsold(lmold+1),nivsigold(lmold))
326	allocate(apold(lmold),bpold(lmold))
327	allocate(presnivsold(lmold))
328	allocate(uold(imold+1,jmold+1,lmold))
329	allocate(vold(imold+1,jmold+1,lmold))
330	allocate(told(imold+1,jmold+1,lmold))
331	allocate(qold(imold+1,jmold+1,lmold,nqtot))
332	allocate(psold(imold+1,jmold+1))
333	allocate(phisold(imold+1,jmold+1))
334	allocate(tsurfold(imold+1,jmold+1))
335	allocate(tsoilold(imold+1,jmold+1,nsoilmx))
336	allocate(albeold(imold+1,jmold+1),radsolold(imold+1,jmold+1))
337	allocate(sollwold(imold+1,jmold+1),solswold(imold+1,jmold+1))
338	allocate(fderold(imold+1,jmold+1))
339	allocate(sollwdownold(imold+1,jmold+1),dlwold(imold+1,jmold+1))
340	allocate(zmeaold(imold+1,jmold+1),zstdold(imold+1,jmold+1))
341	allocate(zsigold(imold+1,jmold+1),zgamold(imold+1,jmold+1))
342	allocate(ztheold(imold+1,jmold+1),zpicold(imold+1,jmold+1))
343	allocate(zvalold(imold+1,jmold+1))
344
345	allocate(var (imold+1,jmold+1,llm))
346	allocate(varp1 (imold+1,jmold+1,llm+1))
347
348	print*,"Initialisations OK"
349
350	c-----------------------------------------------------------------------
351	c Lecture des longitudes et latitudes
352	c-----------------------------------------------------------------------
353	c
354	ierr = NF_INQ_VARID (nid, "rlonv", nvarid)
355	IF (ierr .NE. NF_NOERR) THEN
356	PRINT*, "new_grid: Le champ <rlonv> est absent de start.nc"
357	CALL abort
358	ENDIF
359	#ifdef NC_DOUBLE
360	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlonvold)
361	#else
362	ierr = NF_GET_VAR_REAL(nid, nvarid, rlonvold)
363	#endif
364	IF (ierr .NE. NF_NOERR) THEN
365	PRINT*, "new_grid: Lecture echouee pour <rlonv>"
366	CALL abort
367	ENDIF
368	c
369	ierr = NF_INQ_VARID (nid, "rlatu", nvarid)
370	IF (ierr .NE. NF_NOERR) THEN
371	PRINT*, "new_grid: Le champ <rlatu> est absent de start.nc"
372	CALL abort
373	ENDIF
374	#ifdef NC_DOUBLE
375	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlatuold)
376	#else
377	ierr = NF_GET_VAR_REAL(nid, nvarid, rlatuold)
378	#endif
379	IF (ierr .NE. NF_NOERR) THEN
380	PRINT*, "new_grid: Lecture echouee pour <rlatu>"
381	CALL abort
382	ENDIF
383	c
384	ierr = NF_INQ_VARID (nid, "rlonu", nvarid)
385	IF (ierr .NE. NF_NOERR) THEN
386	PRINT*, "new_grid: Le champ <rlonu> est absent de start.nc"
387	CALL abort
388	ENDIF
389	#ifdef NC_DOUBLE
390	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlonuold)
391	#else
392	ierr = NF_GET_VAR_REAL(nid, nvarid, rlonuold)
393	#endif
394	IF (ierr .NE. NF_NOERR) THEN
395	PRINT*, "new_grid: Lecture echouee pour <rlonu>"
396	CALL abort
397	ENDIF
398	c
399	ierr = NF_INQ_VARID (nid, "rlatv", nvarid)
400	IF (ierr .NE. NF_NOERR) THEN
401	PRINT*, "new_grid: Le champ <rlatv> est absent de start.nc"
402	CALL abort
403	ENDIF
404	#ifdef NC_DOUBLE
405	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlatvold)
406	#else
407	ierr = NF_GET_VAR_REAL(nid, nvarid, rlatvold)
408	#endif
409	IF (ierr .NE. NF_NOERR) THEN
410	PRINT*, "new_grid: Lecture echouee pour <rlatv>"
411	CALL abort
412	ENDIF
413	c
414
415	print*,"Lecture lon/lat OK"
416
417	c-----------------------------------------------------------------------
418	c Lecture des niveaux verticaux
419	c-----------------------------------------------------------------------
420	c
421	ierr = NF_INQ_VARID (nid, "nivsigs", nvarid)
422	IF (ierr .NE. NF_NOERR) THEN
423	PRINT*, "dynetat0: Le champ <nivsigs> est absent"
424	CALL abort
425	ENDIF
426	#ifdef NC_DOUBLE
427	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, nivsigsold)
428	#else
429	ierr = NF_GET_VAR_REAL(nid, nvarid, nivsigsold)
430	#endif
431	IF (ierr .NE. NF_NOERR) THEN
432	PRINT*, "dynetat0: Lecture echouee pour <nivsigs>"
433	CALL abort
434	ENDIF
435
436	ierr = NF_INQ_VARID (nid, "nivsig", nvarid)
437	IF (ierr .NE. NF_NOERR) THEN
438	PRINT*, "dynetat0: Le champ <nivsig> est absent"
439	CALL abort
440	ENDIF
441	#ifdef NC_DOUBLE
442	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, nivsigold)
443	#else
444	ierr = NF_GET_VAR_REAL(nid, nvarid, nivsigold)
445	#endif
446	IF (ierr .NE. NF_NOERR) THEN
447	PRINT*, "dynetat0: Lecture echouee pour <nivsig>"
448	CALL abort
449	ENDIF
450
451	ierr = NF_INQ_VARID (nid, "ap", nvarid)
452	IF (ierr .NE. NF_NOERR) THEN
453	PRINT*, "new_grid: Le champ <ap> est absent de start.nc"
454	CALL abort
455	ELSE
456	#ifdef NC_DOUBLE
457	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, apold)
458	#else
459	ierr = NF_GET_VAR_REAL(nid, nvarid, apold)
460	#endif
461	IF (ierr .NE. NF_NOERR) THEN
462	PRINT*, "new_grid: Lecture echouee pour <ap>"
463	ENDIF
464	ENDIF
465	c
466	ierr = NF_INQ_VARID (nid, "bp", nvarid)
467	IF (ierr .NE. NF_NOERR) THEN
468	PRINT*, "new_grid: Le champ <bp> est absent de start.nc"
469	CALL abort
470	ENDIF
471	#ifdef NC_DOUBLE
472	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, bpold)
473	#else
474	ierr = NF_GET_VAR_REAL(nid, nvarid, bpold)
475	#endif
476	IF (ierr .NE. NF_NOERR) THEN
477	PRINT*, "new_grid: Lecture echouee pour <bp>"
478	CALL abort
479	END IF
480
481	ierr = NF_INQ_VARID (nid, "presnivs", nvarid)
482	IF (ierr .NE. NF_NOERR) THEN
483	PRINT*, "dynetat0: Le champ <presnivs> est absent"
484	CALL abort
485	ENDIF
486	#ifdef NC_DOUBLE
487	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, presnivsold)
488	#else
489	ierr = NF_GET_VAR_REAL(nid, nvarid, presnivsold)
490	#endif
491	IF (ierr .NE. NF_NOERR) THEN
492	PRINT*, "dynetat0: Lecture echouee pour <presnivs>"
493	CALL abort
494	ENDIF
495
496	c-----------------------------------------------------------------------
497	c Lecture geopotentiel au sol
498	c-----------------------------------------------------------------------
499	c
500	ierr = NF_INQ_VARID (nid, "phisinit", nvarid)
501	IF (ierr .NE. NF_NOERR) THEN
502	PRINT*, "new_grid: Le champ <phisinit> est absent de start.nc"
503	CALL abort
504	ENDIF
505	#ifdef NC_DOUBLE
506	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, phisold)
507	#else
508	ierr = NF_GET_VAR_REAL(nid, nvarid, phisold)
509	#endif
510	IF (ierr .NE. NF_NOERR) THEN
511	PRINT*, "new_grid: Lecture echouee pour <phisinit>"
512	CALL abort
513	ENDIF
514
515	print*,"Lecture niveaux et geopot OK"
516
517	c-----------------------------------------------------------------------
518	c Lecture des champs 2D ()
519	c-----------------------------------------------------------------------
520
521	start=(/1,1,1,0/)
522	counter=(/imold+1,jmold+1,1,0/)
523
524	ierr = NF_INQ_VARID (nid, "ps", nvarid)
525	IF (ierr .NE. NF_NOERR) THEN
526	PRINT*, "lect_start_archive: Le champ <ps> est absent"
527	CALL abort
528	ENDIF
529	#ifdef NC_DOUBLE
530	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,psold)
531	#else
532	ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,psold)
533	#endif
534	IF (ierr .NE. NF_NOERR) THEN
535	PRINT*, "lect_start_archive: Lecture echouee pour <ps>"
536	CALL abort
537	ENDIF
538	c
539	ierr = NF_INQ_VARID (nid, "tsurf", nvarid)
540	IF (ierr .NE. NF_NOERR) THEN
541	PRINT*, "lect_start_archive: Le champ <tsurf> est absent"
542	CALL abort
543	ENDIF
544	#ifdef NC_DOUBLE
545	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,tsurfold)
546	#else
547	ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,tsurfold)
548	#endif
549	IF (ierr .NE. NF_NOERR) THEN
550	PRINT*, "lect_start_archive: Lecture echouee pour <tsurf>"
551	CALL abort
552	ENDIF
553	c
554	do isoil=1,nsoilmx
555	write(str2,'(i2.2)') isoil
556	c
557	ierr = NF_INQ_VARID (nid, "Tsoil"//str2, nvarid)
558	IF (ierr .NE. NF_NOERR) THEN
559	PRINT*, "lect_start_archive:
560	. Le champ <","Tsoil"//str2,"> est absent"
561	CALL abort
562	ENDIF
563	#ifdef NC_DOUBLE
564	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,
565	. tsoilold(1,1,isoil))
566	#else
567	ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,
568	. tsoilold(1,1,isoil))
569	#endif
570	IF (ierr .NE. NF_NOERR) THEN
571	PRINT*, "lect_start_archive:
572	. Lecture echouee pour <","Tsoil"//str2,">"
573	CALL abort
574	ENDIF
575	end do
576	c
577	ierr = NF_INQ_VARID (nid, "albe", nvarid)
578	IF (ierr .NE. NF_NOERR) THEN
579	PRINT*, "lect_start_archive: Le champ <albe> est absent"
580	CALL abort
581	ENDIF
582	#ifdef NC_DOUBLE
583	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,albeold)
584	#else
585	ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,albeold)
586	#endif
587	IF (ierr .NE. NF_NOERR) THEN
588	PRINT*, "lect_start_archive: Lecture echouee pour <albe>"
589	CALL abort
590	ENDIF
591	c
592	ierr = NF_INQ_VARID (nid, "radsol", nvarid)
593	IF (ierr .NE. NF_NOERR) THEN
594	PRINT*, "lect_start_archive: Le champ <radsol> est absent"
595	CALL abort
596	ENDIF
597	#ifdef NC_DOUBLE
598	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,radsolold)
599	#else
600	ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,radsolold)
601	#endif
602	IF (ierr .NE. NF_NOERR) THEN
603	PRINT*, "lect_start_archive: Lecture echouee pour <radsol>"
604	CALL abort
605	ENDIF
606	c
607	ierr = NF_INQ_VARID (nid, "sollw", nvarid)
608	IF (ierr .NE. NF_NOERR) THEN
609	PRINT*, "lect_start_archive: Le champ <sollw> est absent"
610	CALL abort
611	ENDIF
612	#ifdef NC_DOUBLE
613	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,sollwold)
614	#else
615	ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,sollwold)
616	#endif
617	IF (ierr .NE. NF_NOERR) THEN
618	PRINT*, "lect_start_archive: Lecture echouee pour <sollw>"
619	CALL abort
620	ENDIF
621	c
622	ierr = NF_INQ_VARID (nid, "solsw", nvarid)
623	IF (ierr .NE. NF_NOERR) THEN
624	PRINT*, "lect_start_archive: Le champ <solsw> est absent"
625	CALL abort
626	ENDIF
627	#ifdef NC_DOUBLE
628	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,solswold)
629	#else
630	ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,solswold)
631	#endif
632	IF (ierr .NE. NF_NOERR) THEN
633	PRINT*, "lect_start_archive: Lecture echouee pour <solsw>"
634	CALL abort
635	ENDIF
636	c
637	ierr = NF_INQ_VARID (nid, "fder", nvarid)
638	IF (ierr .NE. NF_NOERR) THEN
639	PRINT*, "lect_start_archive: Le champ <fder> est absent"
640	CALL abort
641	ENDIF
642	#ifdef NC_DOUBLE
643	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,fderold)
644	#else
645	ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,fderold)
646	#endif
647	IF (ierr .NE. NF_NOERR) THEN
648	PRINT*, "lect_start_archive: Lecture echouee pour <fder>"
649	CALL abort
650	ENDIF
651	c
652	ierr = NF_INQ_VARID (nid, "dlw", nvarid)
653	IF (ierr .NE. NF_NOERR) THEN
654	PRINT*, "lect_start_archive: Le champ <dlw> est absent"
655	CALL abort
656	ENDIF
657	#ifdef NC_DOUBLE
658	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,dlwold)
659	#else
660	ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,dlwold)
661	#endif
662	IF (ierr .NE. NF_NOERR) THEN
663	PRINT*, "lect_start_archive: Lecture echouee pour <dlw>"
664	CALL abort
665	ENDIF
666	c
667	ierr = NF_INQ_VARID (nid, "sollwdown", nvarid)
668	IF (ierr .NE. NF_NOERR) THEN
669	PRINT*, "lect_start_archive: Le champ <sollwdown> est absent"
670	CALL abort
671	ENDIF
672	#ifdef NC_DOUBLE
673	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,sollwdownold)
674	#else
675	ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,sollwdownold)
676	#endif
677	IF (ierr .NE. NF_NOERR) THEN
678	PRINT*, "lect_start_archive: Lecture echouee pour <sollwdown>"
679	CALL abort
680	ENDIF
681	c
682	ierr = NF_INQ_VARID (nid, "zmea", nvarid)
683	IF (ierr .NE. NF_NOERR) THEN
684	PRINT*, "lect_start_archive: Le champ <zmea> est absent"
685	PRINT*, " Il est donc initialise a zero"
686	zmeaold=0.
687	ENDIF
688	#ifdef NC_DOUBLE
689	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,zmeaold)
690	#else
691	ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,zmeaold)
692	#endif
693	IF (ierr .NE. NF_NOERR) THEN
694	PRINT*, "lect_start_archive: Lecture echouee pour <zmea>"
695	CALL abort
696	ENDIF
697	c
698	ierr = NF_INQ_VARID (nid, "zstd", nvarid)
699	IF (ierr .NE. NF_NOERR) THEN
700	PRINT*, "lect_start_archive: Le champ <zstd> est absent"
701	PRINT*, " Il est donc initialise a zero"
702	zstdold=0.
703	ENDIF
704	#ifdef NC_DOUBLE
705	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,zstdold)
706	#else
707	ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,zstdold)
708	#endif
709	IF (ierr .NE. NF_NOERR) THEN
710	PRINT*, "lect_start_archive: Lecture echouee pour <zstd>"
711	CALL abort
712	ENDIF
713	c
714	ierr = NF_INQ_VARID (nid, "zsig", nvarid)
715	IF (ierr .NE. NF_NOERR) THEN
716	PRINT*, "lect_start_archive: Le champ <zsig> est absent"
717	PRINT*, " Il est donc initialise a zero"
718	zsigold=0.
719	ENDIF
720	#ifdef NC_DOUBLE
721	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,zsigold)
722	#else
723	ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,zsigold)
724	#endif
725	IF (ierr .NE. NF_NOERR) THEN
726	PRINT*, "lect_start_archive: Lecture echouee pour <zsig>"
727	CALL abort
728	ENDIF
729	c
730	ierr = NF_INQ_VARID (nid, "zgam", nvarid)
731	IF (ierr .NE. NF_NOERR) THEN
732	PRINT*, "lect_start_archive: Le champ <zgam> est absent"
733	PRINT*, " Il est donc initialise a zero"
734	zgamold=0.
735	ENDIF
736	#ifdef NC_DOUBLE
737	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,zgamold)
738	#else
739	ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,zgamold)
740	#endif
741	IF (ierr .NE. NF_NOERR) THEN
742	PRINT*, "lect_start_archive: Lecture echouee pour <zgam>"
743	CALL abort
744	ENDIF
745	c
746	ierr = NF_INQ_VARID (nid, "zthe", nvarid)
747	IF (ierr .NE. NF_NOERR) THEN
748	PRINT*, "lect_start_archive: Le champ <zthe> est absent"
749	PRINT*, " Il est donc initialise a zero"
750	ztheold=0.
751	ENDIF
752	#ifdef NC_DOUBLE
753	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,ztheold)
754	#else
755	ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,ztheold)
756	#endif
757	IF (ierr .NE. NF_NOERR) THEN
758	PRINT*, "lect_start_archive: Lecture echouee pour <zthe>"
759	CALL abort
760	ENDIF
761	c
762	ierr = NF_INQ_VARID (nid, "zpic", nvarid)
763	IF (ierr .NE. NF_NOERR) THEN
764	PRINT*, "lect_start_archive: Le champ <zpic> est absent"
765	PRINT*, " Il est donc initialise a zero"
766	zpicold=0.
767	ENDIF
768	#ifdef NC_DOUBLE
769	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,zpicold)
770	#else
771	ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,zpicold)
772	#endif
773	IF (ierr .NE. NF_NOERR) THEN
774	PRINT*, "lect_start_archive: Lecture echouee pour <zpic>"
775	CALL abort
776	ENDIF
777	c
778	ierr = NF_INQ_VARID (nid, "zval", nvarid)
779	IF (ierr .NE. NF_NOERR) THEN
780	PRINT*, "lect_start_archive: Le champ <zval> est absent"
781	PRINT*, " Il est donc initialise a zero"
782	zvalold=0.
783	ENDIF
784	#ifdef NC_DOUBLE
785	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,zvalold)
786	#else
787	ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,zvalold)
788	#endif
789	IF (ierr .NE. NF_NOERR) THEN
790	PRINT*, "lect_start_archive: Lecture echouee pour <zval>"
791	CALL abort
792	ENDIF
793	c
794
795	print*,"Lecture champs 2D OK"
796
797	c-----------------------------------------------------------------------
798	c Lecture des champs 3D ()
799	c-----------------------------------------------------------------------
800
801	start=(/1,1,1,1/)
802	counter=(/imold+1,jmold+1,lmold,1/)
803	c
804	ierr = NF_INQ_VARID (nid,"temp", nvarid)
805	IF (ierr .NE. NF_NOERR) THEN
806	PRINT*, "lect_start_archive: Le champ <temp> est absent"
807	CALL abort
808	ENDIF
809	#ifdef NC_DOUBLE
810	ierr = NF_GET_VARA_DOUBLE(nid, nvarid, start, counter, told)
811	#else
812	ierr = NF_GET_VARA_REAL(nid, nvarid, start, counter, told)
813	#endif
814	IF (ierr .NE. NF_NOERR) THEN
815	PRINT*, "lect_start_archive: Lecture echouee pour <temp>"
816	CALL abort
817	ENDIF
818	c
819	ierr = NF_INQ_VARID (nid,"u", nvarid)
820	IF (ierr .NE. NF_NOERR) THEN
821	PRINT*, "lect_start_archive: Le champ <u> est absent"
822	CALL abort
823	ENDIF
824	#ifdef NC_DOUBLE
825	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,uold)
826	#else
827	ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,uold)
828	#endif
829	IF (ierr .NE. NF_NOERR) THEN
830	PRINT*, "lect_start_archive: Lecture echouee pour <u>"
831	CALL abort
832	ENDIF
833	c
834	ierr = NF_INQ_VARID (nid,"v", nvarid)
835	IF (ierr .NE. NF_NOERR) THEN
836	PRINT*, "lect_start_archive: Le champ <v> est absent"
837	CALL abort
838	ENDIF
839	#ifdef NC_DOUBLE
840	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,vold)
841	#else
842	ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,vold)
843	#endif
844	IF (ierr .NE. NF_NOERR) THEN
845	PRINT*, "lect_start_archive: Lecture echouee pour <v>"
846	CALL abort
847	ENDIF
848	c
849
850	c TNAME: IL EST LU A PARTIR DE traceur.def (mettre l'ancien si
851	c changement du nombre de traceurs)
852
853	IF(iflag_trac.eq.1) THEN
854	DO iq=1,nqtot
855	ierr = NF_INQ_VARID (nid, tname(iq), nvarid)
856	IF (ierr .NE. NF_NOERR) THEN
857	PRINT*, "dynetat0: Le champ <"//tname(iq)//"> est absent"
858	PRINT*, " Il est donc initialise a zero"
859	qold(:,:,:,iq)=0.
860	ELSE
861	#ifdef NC_DOUBLE
862	ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,counter,qold(1,1,1,iq))
863	#else
864	ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,qold(1,1,1,iq))
865	#endif
866	IF (ierr .NE. NF_NOERR) THEN
867	PRINT*, "dynetat0: Lecture echouee pour "//tname(iq)
868	CALL abort
869	ENDIF
870	ENDIF
871	ENDDO
872	ENDIF
873
874
875	print*,"Lecture champs 3D OK"
876
877	c=======================================================================
878	c INTERPOLATION DANS LA NOUVELLE GRILLE et initialisation des variables
879	c=======================================================================
880	c Interpolation horizontale puis passage dans la grille physique pour
881	c les variables physique
882	c Interpolation verticale puis horizontale pour chaque variable 3D
883	c=======================================================================
884
885	c-----------------------------------------------------------------------
886	c Variable 2d :
887	c-----------------------------------------------------------------------
888
889	c Relief
890	! topoflag = F: we keep the existing topography
891	! = T: we read it from the Relief.nc file
892	! topoflag need to be in the specific run.def for newstart
893	topoflag = . FALSE .
894	CALL getin('topoflag',topoflag)
895	! notopo = T: we go back to flat surface
896	notopo = .FALSE.
897	CALL getin('notopo',notopo)
898
899	print*,zmeaold(2,1:10)
900
901	IF ( topoflag ) then
902	print*,"Topography (phis) read in file Relief.nc"
903	print*,"For Venus, map directly inverted in Relief.nc"
904	CALL startget_phys2d('surfgeo',iip1,jjp1,rlonv,rlatu,phis,0.0,
905	. jjm ,rlonu,rlatv,.true.)
906	c needed ?
907	phis(iip1,:) = phis(1,:)
908
909	CALL startget_phys1d('zmea',iip1,jjp1,rlonv,rlatu,ngridmx,zmea,
910	. 0.0,jjm,rlonu,rlatv,.true.)
911	CALL startget_phys1d('zstd',iip1,jjp1,rlonv,rlatu,ngridmx,zstd,
912	. 0.0,jjm,rlonu,rlatv,.true.)
913	CALL startget_phys1d('zsig',iip1,jjp1,rlonv,rlatu,ngridmx,zsig,
914	. 0.0,jjm,rlonu,rlatv,.true.)
915	CALL startget_phys1d('zgam',iip1,jjp1,rlonv,rlatu,ngridmx,zgam,
916	. 0.0,jjm,rlonu,rlatv,.true.)
917	CALL startget_phys1d('zthe',iip1,jjp1,rlonv,rlatu,ngridmx,zthe,
918	. 0.0,jjm,rlonu,rlatv,.true.)
919	CALL startget_phys1d('zpic',iip1,jjp1,rlonv,rlatu,ngridmx,zpic,
920	. 0.0,jjm,rlonu,rlatv,.true.)
921	CALL startget_phys1d('zval',iip1,jjp1,rlonv,rlatu,ngridmx,zval,
922	. 0.0,jjm,rlonu,rlatv,.true.)
923
924	ELSE IF ( notopo ) THEN
925	print*,'Flattening the topography'
926	phis=0.
927	zmea=0.
928	zstd=0.
929	zsig=0.
930	zgam=0.
931	zthe=0.
932	zpic=0.
933	zval=0.
934	ELSE
935	print*,'Using existing topography'
936	call interp_horiz (phisold,phis,imold,jmold,iim,jjm,1,
937	& rlonuold,rlatvold,rlonu,rlatv)
938
939	call interp_horiz (zmeaold,zmeaS,imold,jmold,iim,jjm,1,
940	& rlonuold,rlatvold,rlonu,rlatv)
941	call gr_dyn_fi (1,iip1,jjp1,ngridmx,zmeaS,zmea)
942	call interp_horiz (zstdold,zstdS,imold,jmold,iim,jjm,1,
943	& rlonuold,rlatvold,rlonu,rlatv)
944	call gr_dyn_fi (1,iip1,jjp1,ngridmx,zstdS,zstd)
945	call interp_horiz (zsigold,zsigS,imold,jmold,iim,jjm,1,
946	& rlonuold,rlatvold,rlonu,rlatv)
947	call gr_dyn_fi (1,iip1,jjp1,ngridmx,zsigS,zsig)
948	call interp_horiz (zgamold,zgamS,imold,jmold,iim,jjm,1,
949	& rlonuold,rlatvold,rlonu,rlatv)
950	call gr_dyn_fi (1,iip1,jjp1,ngridmx,zgamS,zgam)
951	call interp_horiz (ztheold,ztheS,imold,jmold,iim,jjm,1,
952	& rlonuold,rlatvold,rlonu,rlatv)
953	call gr_dyn_fi (1,iip1,jjp1,ngridmx,ztheS,zthe)
954	call interp_horiz (zpicold,zpicS,imold,jmold,iim,jjm,1,
955	& rlonuold,rlatvold,rlonu,rlatv)
956	call gr_dyn_fi (1,iip1,jjp1,ngridmx,zpicS,zpic)
957	call interp_horiz (zvalold,zvalS,imold,jmold,iim,jjm,1,
958	& rlonuold,rlatvold,rlonu,rlatv)
959	call gr_dyn_fi (1,iip1,jjp1,ngridmx,zvalS,zval)
960	ENDIF
961
962	print*,"New surface geopotential OK"
963
964	c Lat et lon pour physique
965	do i=1,iip1
966	rlatS(i,:)=rlatu(:)*180./pi
967	enddo
968	call gr_dyn_fi (1,iip1,jjp1,ngridmx,rlatS,rlat)
969
970	do j=2,jjm
971	rlonS(:,j)=rlonv(:)*180./pi
972	enddo
973	rlonS(:,1)=0.
974	rlonS(:,jjp1)=0.
975	call gr_dyn_fi (1,iip1,jjp1,ngridmx,rlonS,rlon)
976
977	c Temperature de surface
978	! razTS need to be in the specific run.def for newstart
979	razTS = . FALSE .
980	CALL getin('razTS',razTS)
981
982	if (razTS) then
983	tsurf(:) = 735.
984	else
985	call interp_horiz (tsurfold,tsurfS,imold,jmold,iim,jjm,1,
986	& rlonuold,rlatvold,rlonu,rlatv)
987	call gr_dyn_fi (1,iip1,jjp1,ngridmx,tsurfS,tsurf)
988	c write(44,*) 'tsurf', tsurf
989	endif
990
991	c Temperature du sous-sol
992	if (razTS) then
993	tsoil(:,:)=735.
994	else
995	call interp_horiz(tsoilold,tsoilS,
996	& imold,jmold,iim,jjm,nsoilmx,
997	& rlonuold,rlatvold,rlonu,rlatv)
998	call gr_dyn_fi (nsoilmx,iip1,jjp1,ngridmx,tsoilS,tsoil)
999	c write(45,*) 'tsoil',tsoil
1000	endif
1001
1002	! CHANGING ALBEDO: may be done through run.def
1003	albedoflag = . FALSE .
1004	CALL getin('albedoflag',albedoflag)
1005
1006	IF ( albedoflag ) then
1007	print*,"Albedo is reinitialized to the albedo value in run.def"
1008	print*,"We may want to consider a map later on..."
1009	albedo=0.1
1010	CALL getin('albedo',albedo)
1011	albe=albedo
1012	ELSE
1013	call interp_horiz (albeold,albeS,imold,jmold,iim,jjm,1,
1014	& rlonuold,rlatvold,rlonu,rlatv)
1015	call gr_dyn_fi (1,iip1,jjp1,ngridmx,albeS,albe)
1016	ENDIF
1017
1018	call interp_horiz (radsolold,radsolS,imold,jmold,iim,jjm,1,
1019	& rlonuold,rlatvold,rlonu,rlatv)
1020	call gr_dyn_fi (1,iip1,jjp1,ngridmx,radsolS,radsol)
1021
1022	call interp_horiz (sollwold,sollwS,imold,jmold,iim,jjm,1,
1023	& rlonuold,rlatvold,rlonu,rlatv)
1024	call gr_dyn_fi (1,iip1,jjp1,ngridmx,sollwS,sollw)
1025
1026	call interp_horiz (solswold,solswS,imold,jmold,iim,jjm,1,
1027	& rlonuold,rlatvold,rlonu,rlatv)
1028	call gr_dyn_fi (1,iip1,jjp1,ngridmx,solswS,solsw)
1029
1030	call interp_horiz (fderold,fderS,imold,jmold,iim,jjm,1,
1031	& rlonuold,rlatvold,rlonu,rlatv)
1032	call gr_dyn_fi (1,iip1,jjp1,ngridmx,fderS,fder)
1033
1034	call interp_horiz (dlwold,dlwS,imold,jmold,iim,jjm,1,
1035	& rlonuold,rlatvold,rlonu,rlatv)
1036	call gr_dyn_fi (1,iip1,jjp1,ngridmx,dlwS,dlw)
1037
1038	call interp_horiz (sollwdownold,sollwdownS,imold,jmold,iim,jjm,1,
1039	& rlonuold,rlatvold,rlonu,rlatv)
1040	call gr_dyn_fi (1,iip1,jjp1,ngridmx,sollwdownS,sollwdown)
1041
1042	print*,"Nouvelles var physiques OK"
1043
1044	c-----------------------------------------------------------------------
1045	c Traitement special de la pression au sol :
1046	c-----------------------------------------------------------------------
1047
1048	IF ( notopo ) THEN
1049	ps=9.2e6
1050	ELSE
1051	c Extrapolation la pression dans la nouvelle grille
1052	call interp_horiz(psold,ps,imold,jmold,iim,jjm,1,
1053	& rlonuold,rlatvold,rlonu,rlatv)
1054	ENDIF
1055
1056	c On assure la conservation de la masse de l'atmosphere
1057	c--------------------------------------------------------------
1058
1059	ptotal = 0.
1060	DO j=1,jjp1
1061	DO i=1,iim
1062	ptotal=ptotal+ps(i,j)*aire(i,j)/g
1063	END DO
1064	END DO
1065
1066	write(,)
1067	write(,)'Ancienne grille: masse de l atm :',ptotalold
1068	write(,)'Nouvelle grille: masse de l atm :',ptotal
1069	write (,) 'Ratio new atm./ old atm =', ptotal/ptotalold
1070	write(,)
1071
1072
1073	DO j=1,jjp1
1074	DO i=1,iip1
1075	ps(i,j)=ps(i,j) * ptotalold/ptotal
1076	END DO
1077	END DO
1078
1079	c la pression de surface et les temperatures ne sont pas reequilibrees en fonction
1080	c de la nouvelle topographie...
1081	c Si l'ajustement inevitable du debut pose des problemes, voir le newstart martien.
1082
1083	print*,"Nouvelle ps OK"
1084	print*,"If changes done on topography, beware !"
1085	print*,"Some time may be needed for adjustments at the beginning"
1086	print*,"so if unstable, relax the timestep and/or dissipation."
1087
1088	c-----------------------------------------------------------------------
1089	c Variable 3d :
1090	c-----------------------------------------------------------------------
1091
1092	CALL pression(ip1jmp1, ap, bp, ps, p3d)
1093	if (disvert_type==1) then
1094	CALL exner_hyb( ip1jmp1, ps, p3d,pks, pk, pkf )
1095	else ! we assume that we are in the disvert_type==2 case
1096	CALL exner_milieu( ip1jmp1, ps, p3d, pks, pk, pkf )
1097	endif
1098
1099	c temperatures atmospheriques
1100
1101	c ATTENTION: peut servir, mais bon...
1102	c modif: profil uniforme
1103	c do l=1,lmold
1104	c do j=1,jmold+1
1105	c do i=1,imold+1
1106	c told(i,j,l)=told(1,jmold/2,l)
1107	c enddo
1108	c enddo
1109	c enddo
1110
1111	! raztemp need to be in the specific run.def for newstart
1112	raztemp = . FALSE .
1113	CALL getin('raztemp',raztemp)
1114
1115	! Reinitialisation of temperature to VIRA profile lisse
1116	if (raztemp) then
1117
1118	allocate(tvira(0:lmold),dzst(0:lmold),zkm(0:lmold))
1119	print*,"Venus = temperature initiale imposee = VIRA lisse "
1120	dzst(0) = 0.0
1121	dzst(1) = -log(p3d(1,1,2)/preff)*r/g
1122	do l=2,lmold
1123	dzst(l)=-log(p3d(1,1,l+1)/p3d(1,1,l))*r/g
1124	enddo
1125	tvira(0) = 735.
1126	zkm(0) = 0.0
1127	do l=1,lmold
1128	zkm(l) = zkm(l-1)+tvira(l-1)*dzst(l)/1000. ! approx avec T(l-1)
1129	if(zkm(l).lt.60.) then
1130	tvira(l)=735.-7.95*zkm(l)
1131	else
1132	tvira(l)=AMAX1(258.-3.*(zkm(l)-60.),168.)
1133	endif
1134	zkm(l) = zkm(l-1)+(tvira(l-1)+tvira(l))/2.*dzst(l)/1000.
1135	enddo
1136	do l=1,lmold
1137	do j=1,jmold+1
1138	do i=1,imold+1
1139	told(i,j,l)=tvira(l)
1140	enddo
1141	enddo
1142	enddo
1143	endif ! end raztemp
1144
1145	write (,) 'told ', told (1,jmold+1,1) ! INFO
1146	call interp_vert
1147	& (told,var,lmold,llm,apold,bpold,ap,bp,
1148	& psold,(imold+1)*(jmold+1))
1149	write (,) 'var ', var (1,jmold+1,1) ! INFO
1150	call interp_horiz(var,t,imold,jmold,iim,jjm,llm,
1151	& rlonuold,rlatvold,rlonu,rlatv)
1152	write (,) 'T ', T(1,jjp1,1) ! INFO
1153	! pour info:
1154	! Si extension verticale, la T est extrapolee constante au-dessus de lmold
1155
1156	c passage grille physique pour restartphy.nc
1157	call gr_dyn_fi (llm,iip1,jjp1,ngridmx,T,t_fi)
1158
1159	! ADAPTATION GCM POUR CP(T)
1160	c passage en temperature potentielle
1161	call t2tpot(ip1jmp1*llm,T,teta,pk)
1162	c on assure la periodicite
1163	teta(iip1,:,:) = teta(1,:,:)
1164
1165	! RESETING U TO uini: may be done through run.def
1166	uini = .FALSE.
1167	CALL getin('uini',uini)
1168	! RESETING U TO 0: may be done through run.def
1169	razvitu = .FALSE.
1170	CALL getin('razvitu',razvitu)
1171	razvitv = .FALSE.
1172	CALL getin('razvitv',razvitv)
1173
1174	c calcul des champ de vent; passage en vent covariant
1175	write (,) 'uold ', uold (1,2,1) ! INFO
1176	call interp_vert
1177	& (uold,var,lmold,llm,apold,bpold,ap,bp,
1178	& psold,(imold+1)*(jmold+1))
1179	write (,) 'var ', var (1,2,1) ! INFO
1180	call interp_horiz(var,us,imold,jmold,iim,jjm,llm,
1181	& rlonuold,rlatvold,rlonu,rlatv)
1182	write (,) 'us ', us (1,2,1) ! INFO
1183
1184	call interp_vert
1185	& (vold,var,lmold,llm,
1186	& apold,bpold,ap,bp,psold,(imold+1)*(jmold+1))
1187	call interp_horiz(var,vs,imold,jmold,iim,jjm,llm,
1188	& rlonuold,rlatvold,rlonu,rlatv)
1189	call scal_wind(us,vs,unat,vnat)
1190	! Reseting u=0
1191	if ((razvitu).and..not.(uini)) then
1192	unat(:,:,:) = 0.
1193	endif
1194	! Reseting u=uini
1195	if ((uini).and..not.(razvitu)) then
1196	do j=1,jjp1
1197	do l=1,llm
1198	if (p3d(1,j,l).gt.3e3) then
1199	unat(:,j,l) = -110./8.03*log(p3d(:,j,l)/9.2e6)
1200	else
1201	unat(:,j,l) = 110./6.62*(log(p3d(:,j,l)/9.2e6)+14.65)
1202	endif
1203	if (abs(rlatS(1,j)).gt.50.) then
1204	unat(:,j,l)=unat(:,j,l)*(90.-abs(rlatS(:,j)))/40.
1205	endif
1206	enddo
1207	enddo
1208	endif
1209	! incompatible options
1210	if ((uini).and.(razvitu)) then
1211	print*,"You have to choose between razvitu and uini..."
1212	stop
1213	endif
1214	write (,) 'unat ', unat (1,2,1) ! INFO
1215
1216	do l=1,llm
1217	do j = 1, jjp1
1218	do i=1,iip1
1219	ucov( i,j,l ) = unat( i,j,l ) * cu(i,j)
1220	! pour info:
1221	! Si extension verticale, on impose u=0 au-dessus de lmold
1222	! if (presnivs(l).lt.presnivsold(lmold)) ucov( i,j,l ) = 0
1223	end do
1224	end do
1225	end do
1226	write (,) 'ucov ', ucov (1,2,1) ! INFO
1227	c write(48,*) 'ucov',ucov
1228	! Reseting v=0
1229	if (razvitv) then
1230	vnat(:,:,:) = 0.
1231	endif
1232	write (,) 'vnat ', vnat (1,2,1) ! INFO
1233	do l=1,llm
1234	do j = 1, jjm
1235	do i=1,iim
1236	vcov( i,j,l ) = vnat( i,j,l ) * cv(i,j)
1237	! pour info:
1238	! Si extension verticale, on impose v=0 au-dessus de lmold
1239	! if (presnivs(l).lt.presnivsold(lmold)) vcov( i,j,l ) = 0
1240	end do
1241	vcov( iip1,j,l ) = vcov( 1,j,l )
1242	end do
1243	end do
1244	c write(49,*) 'ucov',vcov
1245
1246	c masse: on la recalcule (ps a été ajustée pour conserver la masse totale)
1247	call massdair(p3d,masse)
1248
1249	c traceurs 3D
1250	do iq = 1, nqtot
1251	call interp_vert(qold(1,1,1,iq),var,lmold,llm,
1252	& apold,bpold,ap,bp,psold,(imold+1)*(jmold+1))
1253	call interp_horiz(var,q(1,1,1,iq),imold,jmold,iim,jjm,llm,
1254	& rlonuold,rlatvold,rlonu,rlatv)
1255	enddo
1256
1257	print*,"Nouvelles var dynamiques OK"
1258
1259	c=======================================================================
1260	c Ecriture des restart.nc et restartphy.nc :
1261	c=======================================================================
1262
1263	call writerestart('restart.nc',tab_cntrl_dyn,
1264	. phis,vcov,ucov,teta,masse,q,ps)
1265
1266	print*,"restart OK"
1267
1268	call writerestartphy('restartphy.nc',tab_cntrl_fi,ngridmx,llm,
1269	. rlat,rlon,tsurf,tsoil,albe,solsw, sollw,fder,dlw,
1270	. sollwdown,radsol,
1271	. zmea,zstd,zsig,zgam,zthe,zpic,zval,
1272	. t_fi)
1273
1274	print*,"restartphy OK"
1275
1276	end

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

Download in other formats:

Original Format