Context Navigation

newstart.F @ 2013

Last change on this file since 2013 was 1523, checked in by emillour, 9 years ago

All models: More updates to make planetary codes (+Earth) setups converge.

in dyn3d_common:
convmas.F => convmas.F90
enercin.F => enercin.F90
flumass.F => flumass.F90
massbar.F => massbar.F90
tourpot.F => tourpot.F90
vitvert.F => vitvert.F90

in misc:
move "q_sat" from "dyn3d_common" to "misc" (in Earth model, it is also called by the physics)
move "write_field" from "dyn3d_common" to "misc"(may be called from physics or dynamics and depends on neither).

in phy_common:
move "write_field_phy" here since it may be called from any physics package)
add module "regular_lonlat_mod" to store global information on lon-lat grid

in dynlonlat_phylonlat/phy*:
turn "iniphysiq.F90" into module "iniphysiq_mod.F90" (and of course adapt gcm.F[90] and 1D models accordingly)

File size: 42.0 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	call ini_cpdet
253
254	c=======================================================================
255	c CHANGEMENT DE CONSTANTES CONTENUES DANS tab_cntrl
256	c=======================================================================
257	c changement de la resolution dans le fichier de controle
258	tab_cntrl(1) = REAL(iim)
259	tab_cntrl(2) = REAL(jjm)
260	tab_cntrl(3) = REAL(llm)
261
262	c--------------------------------
263	c We use a specific run.def to read new parameters that need to be changed
264	c--------------------------------
265
266	c Changement de cpp:
267	CALL getin('cpp',cpp)
268	kappa = (8.314511/43.44E-3)/cpp
269	tab_cntrl(9) = cpp
270	tab_cntrl(10) = kappa
271
272	c CHANGING THE ZOOM PARAMETERS TO CHANGE THE GRID
273	CALL getin('clon',clon)
274	CALL getin('clat',clat)
275	tab_cntrl(20) = clon
276	tab_cntrl(21) = clat
277	CALL getin('grossismx',grossismx)
278	CALL getin('grossismy',grossismy)
279	tab_cntrl(22) = grossismx
280	tab_cntrl(23) = grossismy
281	CALL getin('fxyhypb',fxyhypb)
282	IF ( fxyhypb ) THEN
283	CALL getin('dzoomx',dzoomx)
284	CALL getin('dzoomy',dzoomy)
285	tab_cntrl(25) = dzoomx
286	tab_cntrl(26) = dzoomy
287	CALL getin('taux',taux)
288	CALL getin('tauy',tauy)
289	tab_cntrl(28) = taux
290	tab_cntrl(29) = tauy
291	ELSE
292	CALL getin('ysinus',ysinus)
293	IF (ysinus) THEN
294	tab_cntrl(27) = 1
295	ELSE
296	tab_cntrl(27) = 0
297	ENDIF
298	ENDIF
299
300	c changes must be done BEFORE these lines...
301	DO l=1,length
302	tab_cntrl_dyn(l)= tab_cntrl(l)
303	tab_cntrl_fi(l) = tab_cntrl(l+length)
304	ENDDO
305
306	c-----------------------------------------------------------------------
307	c Autres initialisations
308	c-----------------------------------------------------------------------
309
310	CALL iniconst
311	CALL inigeom
312	call inifilr
313
314	c-----------------------------------------------------------------------
315	c Allocations des anciennes variables
316	c-----------------------------------------------------------------------
317
318	allocate(rlonuold(imold+1), rlatvold(jmold))
319	allocate(rlonvold(imold+1), rlatuold(jmold+1))
320	allocate(nivsigsold(lmold+1),nivsigold(lmold))
321	allocate(apold(lmold),bpold(lmold))
322	allocate(presnivsold(lmold))
323	allocate(uold(imold+1,jmold+1,lmold))
324	allocate(vold(imold+1,jmold+1,lmold))
325	allocate(told(imold+1,jmold+1,lmold))
326	allocate(qold(imold+1,jmold+1,lmold,nqtot))
327	allocate(psold(imold+1,jmold+1))
328	allocate(phisold(imold+1,jmold+1))
329	allocate(tsurfold(imold+1,jmold+1))
330	allocate(tsoilold(imold+1,jmold+1,nsoilmx))
331	allocate(albeold(imold+1,jmold+1),radsolold(imold+1,jmold+1))
332	allocate(sollwold(imold+1,jmold+1),solswold(imold+1,jmold+1))
333	allocate(fderold(imold+1,jmold+1))
334	allocate(sollwdownold(imold+1,jmold+1),dlwold(imold+1,jmold+1))
335	allocate(zmeaold(imold+1,jmold+1),zstdold(imold+1,jmold+1))
336	allocate(zsigold(imold+1,jmold+1),zgamold(imold+1,jmold+1))
337	allocate(ztheold(imold+1,jmold+1),zpicold(imold+1,jmold+1))
338	allocate(zvalold(imold+1,jmold+1))
339
340	allocate(var (imold+1,jmold+1,llm))
341	allocate(varp1 (imold+1,jmold+1,llm+1))
342
343	print*,"Initialisations OK"
344
345	c-----------------------------------------------------------------------
346	c Lecture des longitudes et latitudes
347	c-----------------------------------------------------------------------
348	c
349	ierr = NF_INQ_VARID (nid, "rlonv", nvarid)
350	IF (ierr .NE. NF_NOERR) THEN
351	PRINT*, "new_grid: Le champ <rlonv> est absent de start.nc"
352	CALL abort
353	ENDIF
354	#ifdef NC_DOUBLE
355	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlonvold)
356	#else
357	ierr = NF_GET_VAR_REAL(nid, nvarid, rlonvold)
358	#endif
359	IF (ierr .NE. NF_NOERR) THEN
360	PRINT*, "new_grid: Lecture echouee pour <rlonv>"
361	CALL abort
362	ENDIF
363	c
364	ierr = NF_INQ_VARID (nid, "rlatu", nvarid)
365	IF (ierr .NE. NF_NOERR) THEN
366	PRINT*, "new_grid: Le champ <rlatu> est absent de start.nc"
367	CALL abort
368	ENDIF
369	#ifdef NC_DOUBLE
370	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlatuold)
371	#else
372	ierr = NF_GET_VAR_REAL(nid, nvarid, rlatuold)
373	#endif
374	IF (ierr .NE. NF_NOERR) THEN
375	PRINT*, "new_grid: Lecture echouee pour <rlatu>"
376	CALL abort
377	ENDIF
378	c
379	ierr = NF_INQ_VARID (nid, "rlonu", nvarid)
380	IF (ierr .NE. NF_NOERR) THEN
381	PRINT*, "new_grid: Le champ <rlonu> est absent de start.nc"
382	CALL abort
383	ENDIF
384	#ifdef NC_DOUBLE
385	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlonuold)
386	#else
387	ierr = NF_GET_VAR_REAL(nid, nvarid, rlonuold)
388	#endif
389	IF (ierr .NE. NF_NOERR) THEN
390	PRINT*, "new_grid: Lecture echouee pour <rlonu>"
391	CALL abort
392	ENDIF
393	c
394	ierr = NF_INQ_VARID (nid, "rlatv", nvarid)
395	IF (ierr .NE. NF_NOERR) THEN
396	PRINT*, "new_grid: Le champ <rlatv> est absent de start.nc"
397	CALL abort
398	ENDIF
399	#ifdef NC_DOUBLE
400	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlatvold)
401	#else
402	ierr = NF_GET_VAR_REAL(nid, nvarid, rlatvold)
403	#endif
404	IF (ierr .NE. NF_NOERR) THEN
405	PRINT*, "new_grid: Lecture echouee pour <rlatv>"
406	CALL abort
407	ENDIF
408	c
409
410	print*,"Lecture lon/lat OK"
411
412	c-----------------------------------------------------------------------
413	c Lecture des niveaux verticaux
414	c-----------------------------------------------------------------------
415	c
416	ierr = NF_INQ_VARID (nid, "nivsigs", nvarid)
417	IF (ierr .NE. NF_NOERR) THEN
418	PRINT*, "dynetat0: Le champ <nivsigs> est absent"
419	CALL abort
420	ENDIF
421	#ifdef NC_DOUBLE
422	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, nivsigsold)
423	#else
424	ierr = NF_GET_VAR_REAL(nid, nvarid, nivsigsold)
425	#endif
426	IF (ierr .NE. NF_NOERR) THEN
427	PRINT*, "dynetat0: Lecture echouee pour <nivsigs>"
428	CALL abort
429	ENDIF
430
431	ierr = NF_INQ_VARID (nid, "nivsig", nvarid)
432	IF (ierr .NE. NF_NOERR) THEN
433	PRINT*, "dynetat0: Le champ <nivsig> est absent"
434	CALL abort
435	ENDIF
436	#ifdef NC_DOUBLE
437	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, nivsigold)
438	#else
439	ierr = NF_GET_VAR_REAL(nid, nvarid, nivsigold)
440	#endif
441	IF (ierr .NE. NF_NOERR) THEN
442	PRINT*, "dynetat0: Lecture echouee pour <nivsig>"
443	CALL abort
444	ENDIF
445
446	ierr = NF_INQ_VARID (nid, "ap", nvarid)
447	IF (ierr .NE. NF_NOERR) THEN
448	PRINT*, "new_grid: Le champ <ap> est absent de start.nc"
449	CALL abort
450	ELSE
451	#ifdef NC_DOUBLE
452	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, apold)
453	#else
454	ierr = NF_GET_VAR_REAL(nid, nvarid, apold)
455	#endif
456	IF (ierr .NE. NF_NOERR) THEN
457	PRINT*, "new_grid: Lecture echouee pour <ap>"
458	ENDIF
459	ENDIF
460	c
461	ierr = NF_INQ_VARID (nid, "bp", nvarid)
462	IF (ierr .NE. NF_NOERR) THEN
463	PRINT*, "new_grid: Le champ <bp> est absent de start.nc"
464	CALL abort
465	ENDIF
466	#ifdef NC_DOUBLE
467	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, bpold)
468	#else
469	ierr = NF_GET_VAR_REAL(nid, nvarid, bpold)
470	#endif
471	IF (ierr .NE. NF_NOERR) THEN
472	PRINT*, "new_grid: Lecture echouee pour <bp>"
473	CALL abort
474	END IF
475
476	ierr = NF_INQ_VARID (nid, "presnivs", nvarid)
477	IF (ierr .NE. NF_NOERR) THEN
478	PRINT*, "dynetat0: Le champ <presnivs> est absent"
479	CALL abort
480	ENDIF
481	#ifdef NC_DOUBLE
482	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, presnivsold)
483	#else
484	ierr = NF_GET_VAR_REAL(nid, nvarid, presnivsold)
485	#endif
486	IF (ierr .NE. NF_NOERR) THEN
487	PRINT*, "dynetat0: Lecture echouee pour <presnivs>"
488	CALL abort
489	ENDIF
490
491	c-----------------------------------------------------------------------
492	c Lecture geopotentiel au sol
493	c-----------------------------------------------------------------------
494	c
495	ierr = NF_INQ_VARID (nid, "phisinit", nvarid)
496	IF (ierr .NE. NF_NOERR) THEN
497	PRINT*, "new_grid: Le champ <phisinit> est absent de start.nc"
498	CALL abort
499	ENDIF
500	#ifdef NC_DOUBLE
501	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, phisold)
502	#else
503	ierr = NF_GET_VAR_REAL(nid, nvarid, phisold)
504	#endif
505	IF (ierr .NE. NF_NOERR) THEN
506	PRINT*, "new_grid: Lecture echouee pour <phisinit>"
507	CALL abort
508	ENDIF
509
510	print*,"Lecture niveaux et geopot OK"
511
512	c-----------------------------------------------------------------------
513	c Lecture des champs 2D ()
514	c-----------------------------------------------------------------------
515
516	start=(/1,1,1,0/)
517	counter=(/imold+1,jmold+1,1,0/)
518
519	ierr = NF_INQ_VARID (nid, "ps", nvarid)
520	IF (ierr .NE. NF_NOERR) THEN
521	PRINT*, "lect_start_archive: Le champ <ps> est absent"
522	CALL abort
523	ENDIF
524	#ifdef NC_DOUBLE
525	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,psold)
526	#else
527	ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,psold)
528	#endif
529	IF (ierr .NE. NF_NOERR) THEN
530	PRINT*, "lect_start_archive: Lecture echouee pour <ps>"
531	CALL abort
532	ENDIF
533	c
534	ierr = NF_INQ_VARID (nid, "tsurf", nvarid)
535	IF (ierr .NE. NF_NOERR) THEN
536	PRINT*, "lect_start_archive: Le champ <tsurf> est absent"
537	CALL abort
538	ENDIF
539	#ifdef NC_DOUBLE
540	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,tsurfold)
541	#else
542	ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,tsurfold)
543	#endif
544	IF (ierr .NE. NF_NOERR) THEN
545	PRINT*, "lect_start_archive: Lecture echouee pour <tsurf>"
546	CALL abort
547	ENDIF
548	c
549	do isoil=1,nsoilmx
550	write(str2,'(i2.2)') isoil
551	c
552	ierr = NF_INQ_VARID (nid, "Tsoil"//str2, nvarid)
553	IF (ierr .NE. NF_NOERR) THEN
554	PRINT*, "lect_start_archive:
555	. Le champ <","Tsoil"//str2,"> est absent"
556	CALL abort
557	ENDIF
558	#ifdef NC_DOUBLE
559	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,
560	. tsoilold(1,1,isoil))
561	#else
562	ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,
563	. tsoilold(1,1,isoil))
564	#endif
565	IF (ierr .NE. NF_NOERR) THEN
566	PRINT*, "lect_start_archive:
567	. Lecture echouee pour <","Tsoil"//str2,">"
568	CALL abort
569	ENDIF
570	end do
571	c
572	ierr = NF_INQ_VARID (nid, "albe", nvarid)
573	IF (ierr .NE. NF_NOERR) THEN
574	PRINT*, "lect_start_archive: Le champ <albe> est absent"
575	CALL abort
576	ENDIF
577	#ifdef NC_DOUBLE
578	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,albeold)
579	#else
580	ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,albeold)
581	#endif
582	IF (ierr .NE. NF_NOERR) THEN
583	PRINT*, "lect_start_archive: Lecture echouee pour <albe>"
584	CALL abort
585	ENDIF
586	c
587	ierr = NF_INQ_VARID (nid, "radsol", nvarid)
588	IF (ierr .NE. NF_NOERR) THEN
589	PRINT*, "lect_start_archive: Le champ <radsol> est absent"
590	CALL abort
591	ENDIF
592	#ifdef NC_DOUBLE
593	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,radsolold)
594	#else
595	ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,radsolold)
596	#endif
597	IF (ierr .NE. NF_NOERR) THEN
598	PRINT*, "lect_start_archive: Lecture echouee pour <radsol>"
599	CALL abort
600	ENDIF
601	c
602	ierr = NF_INQ_VARID (nid, "sollw", nvarid)
603	IF (ierr .NE. NF_NOERR) THEN
604	PRINT*, "lect_start_archive: Le champ <sollw> est absent"
605	CALL abort
606	ENDIF
607	#ifdef NC_DOUBLE
608	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,sollwold)
609	#else
610	ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,sollwold)
611	#endif
612	IF (ierr .NE. NF_NOERR) THEN
613	PRINT*, "lect_start_archive: Lecture echouee pour <sollw>"
614	CALL abort
615	ENDIF
616	c
617	ierr = NF_INQ_VARID (nid, "solsw", nvarid)
618	IF (ierr .NE. NF_NOERR) THEN
619	PRINT*, "lect_start_archive: Le champ <solsw> est absent"
620	CALL abort
621	ENDIF
622	#ifdef NC_DOUBLE
623	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,solswold)
624	#else
625	ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,solswold)
626	#endif
627	IF (ierr .NE. NF_NOERR) THEN
628	PRINT*, "lect_start_archive: Lecture echouee pour <solsw>"
629	CALL abort
630	ENDIF
631	c
632	ierr = NF_INQ_VARID (nid, "fder", nvarid)
633	IF (ierr .NE. NF_NOERR) THEN
634	PRINT*, "lect_start_archive: Le champ <fder> est absent"
635	CALL abort
636	ENDIF
637	#ifdef NC_DOUBLE
638	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,fderold)
639	#else
640	ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,fderold)
641	#endif
642	IF (ierr .NE. NF_NOERR) THEN
643	PRINT*, "lect_start_archive: Lecture echouee pour <fder>"
644	CALL abort
645	ENDIF
646	c
647	ierr = NF_INQ_VARID (nid, "dlw", nvarid)
648	IF (ierr .NE. NF_NOERR) THEN
649	PRINT*, "lect_start_archive: Le champ <dlw> est absent"
650	CALL abort
651	ENDIF
652	#ifdef NC_DOUBLE
653	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,dlwold)
654	#else
655	ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,dlwold)
656	#endif
657	IF (ierr .NE. NF_NOERR) THEN
658	PRINT*, "lect_start_archive: Lecture echouee pour <dlw>"
659	CALL abort
660	ENDIF
661	c
662	ierr = NF_INQ_VARID (nid, "sollwdown", nvarid)
663	IF (ierr .NE. NF_NOERR) THEN
664	PRINT*, "lect_start_archive: Le champ <sollwdown> est absent"
665	CALL abort
666	ENDIF
667	#ifdef NC_DOUBLE
668	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,sollwdownold)
669	#else
670	ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,sollwdownold)
671	#endif
672	IF (ierr .NE. NF_NOERR) THEN
673	PRINT*, "lect_start_archive: Lecture echouee pour <sollwdown>"
674	CALL abort
675	ENDIF
676	c
677	ierr = NF_INQ_VARID (nid, "zmea", nvarid)
678	IF (ierr .NE. NF_NOERR) THEN
679	PRINT*, "lect_start_archive: Le champ <zmea> est absent"
680	PRINT*, " Il est donc initialise a zero"
681	zmeaold=0.
682	ENDIF
683	#ifdef NC_DOUBLE
684	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,zmeaold)
685	#else
686	ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,zmeaold)
687	#endif
688	IF (ierr .NE. NF_NOERR) THEN
689	PRINT*, "lect_start_archive: Lecture echouee pour <zmea>"
690	CALL abort
691	ENDIF
692	c
693	ierr = NF_INQ_VARID (nid, "zstd", nvarid)
694	IF (ierr .NE. NF_NOERR) THEN
695	PRINT*, "lect_start_archive: Le champ <zstd> est absent"
696	PRINT*, " Il est donc initialise a zero"
697	zstdold=0.
698	ENDIF
699	#ifdef NC_DOUBLE
700	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,zstdold)
701	#else
702	ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,zstdold)
703	#endif
704	IF (ierr .NE. NF_NOERR) THEN
705	PRINT*, "lect_start_archive: Lecture echouee pour <zstd>"
706	CALL abort
707	ENDIF
708	c
709	ierr = NF_INQ_VARID (nid, "zsig", nvarid)
710	IF (ierr .NE. NF_NOERR) THEN
711	PRINT*, "lect_start_archive: Le champ <zsig> est absent"
712	PRINT*, " Il est donc initialise a zero"
713	zsigold=0.
714	ENDIF
715	#ifdef NC_DOUBLE
716	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,zsigold)
717	#else
718	ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,zsigold)
719	#endif
720	IF (ierr .NE. NF_NOERR) THEN
721	PRINT*, "lect_start_archive: Lecture echouee pour <zsig>"
722	CALL abort
723	ENDIF
724	c
725	ierr = NF_INQ_VARID (nid, "zgam", nvarid)
726	IF (ierr .NE. NF_NOERR) THEN
727	PRINT*, "lect_start_archive: Le champ <zgam> est absent"
728	PRINT*, " Il est donc initialise a zero"
729	zgamold=0.
730	ENDIF
731	#ifdef NC_DOUBLE
732	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,zgamold)
733	#else
734	ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,zgamold)
735	#endif
736	IF (ierr .NE. NF_NOERR) THEN
737	PRINT*, "lect_start_archive: Lecture echouee pour <zgam>"
738	CALL abort
739	ENDIF
740	c
741	ierr = NF_INQ_VARID (nid, "zthe", nvarid)
742	IF (ierr .NE. NF_NOERR) THEN
743	PRINT*, "lect_start_archive: Le champ <zthe> est absent"
744	PRINT*, " Il est donc initialise a zero"
745	ztheold=0.
746	ENDIF
747	#ifdef NC_DOUBLE
748	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,ztheold)
749	#else
750	ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,ztheold)
751	#endif
752	IF (ierr .NE. NF_NOERR) THEN
753	PRINT*, "lect_start_archive: Lecture echouee pour <zthe>"
754	CALL abort
755	ENDIF
756	c
757	ierr = NF_INQ_VARID (nid, "zpic", nvarid)
758	IF (ierr .NE. NF_NOERR) THEN
759	PRINT*, "lect_start_archive: Le champ <zpic> est absent"
760	PRINT*, " Il est donc initialise a zero"
761	zpicold=0.
762	ENDIF
763	#ifdef NC_DOUBLE
764	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,zpicold)
765	#else
766	ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,zpicold)
767	#endif
768	IF (ierr .NE. NF_NOERR) THEN
769	PRINT*, "lect_start_archive: Lecture echouee pour <zpic>"
770	CALL abort
771	ENDIF
772	c
773	ierr = NF_INQ_VARID (nid, "zval", nvarid)
774	IF (ierr .NE. NF_NOERR) THEN
775	PRINT*, "lect_start_archive: Le champ <zval> est absent"
776	PRINT*, " Il est donc initialise a zero"
777	zvalold=0.
778	ENDIF
779	#ifdef NC_DOUBLE
780	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,zvalold)
781	#else
782	ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,zvalold)
783	#endif
784	IF (ierr .NE. NF_NOERR) THEN
785	PRINT*, "lect_start_archive: Lecture echouee pour <zval>"
786	CALL abort
787	ENDIF
788	c
789
790	print*,"Lecture champs 2D OK"
791
792	c-----------------------------------------------------------------------
793	c Lecture des champs 3D ()
794	c-----------------------------------------------------------------------
795
796	start=(/1,1,1,1/)
797	counter=(/imold+1,jmold+1,lmold,1/)
798	c
799	ierr = NF_INQ_VARID (nid,"temp", nvarid)
800	IF (ierr .NE. NF_NOERR) THEN
801	PRINT*, "lect_start_archive: Le champ <temp> est absent"
802	CALL abort
803	ENDIF
804	#ifdef NC_DOUBLE
805	ierr = NF_GET_VARA_DOUBLE(nid, nvarid, start, counter, told)
806	#else
807	ierr = NF_GET_VARA_REAL(nid, nvarid, start, counter, told)
808	#endif
809	IF (ierr .NE. NF_NOERR) THEN
810	PRINT*, "lect_start_archive: Lecture echouee pour <temp>"
811	CALL abort
812	ENDIF
813	c
814	ierr = NF_INQ_VARID (nid,"u", nvarid)
815	IF (ierr .NE. NF_NOERR) THEN
816	PRINT*, "lect_start_archive: Le champ <u> est absent"
817	CALL abort
818	ENDIF
819	#ifdef NC_DOUBLE
820	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,uold)
821	#else
822	ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,uold)
823	#endif
824	IF (ierr .NE. NF_NOERR) THEN
825	PRINT*, "lect_start_archive: Lecture echouee pour <u>"
826	CALL abort
827	ENDIF
828	c
829	ierr = NF_INQ_VARID (nid,"v", nvarid)
830	IF (ierr .NE. NF_NOERR) THEN
831	PRINT*, "lect_start_archive: Le champ <v> est absent"
832	CALL abort
833	ENDIF
834	#ifdef NC_DOUBLE
835	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,vold)
836	#else
837	ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,vold)
838	#endif
839	IF (ierr .NE. NF_NOERR) THEN
840	PRINT*, "lect_start_archive: Lecture echouee pour <v>"
841	CALL abort
842	ENDIF
843	c
844
845	c TNAME: IL EST LU A PARTIR DE traceur.def (mettre l'ancien si
846	c changement du nombre de traceurs)
847
848	IF(iflag_trac.eq.1) THEN
849	DO iq=1,nqtot
850	ierr = NF_INQ_VARID (nid, tname(iq), nvarid)
851	IF (ierr .NE. NF_NOERR) THEN
852	PRINT*, "dynetat0: Le champ <"//tname(iq)//"> est absent"
853	PRINT*, " Il est donc initialise a zero"
854	qold(:,:,:,iq)=0.
855	ELSE
856	#ifdef NC_DOUBLE
857	ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,counter,qold(1,1,1,iq))
858	#else
859	ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,qold(1,1,1,iq))
860	#endif
861	IF (ierr .NE. NF_NOERR) THEN
862	PRINT*, "dynetat0: Lecture echouee pour "//tname(iq)
863	CALL abort
864	ENDIF
865	ENDIF
866	ENDDO
867	ENDIF
868
869
870	print*,"Lecture champs 3D OK"
871
872	c=======================================================================
873	c INTERPOLATION DANS LA NOUVELLE GRILLE et initialisation des variables
874	c=======================================================================
875	c Interpolation horizontale puis passage dans la grille physique pour
876	c les variables physique
877	c Interpolation verticale puis horizontale pour chaque variable 3D
878	c=======================================================================
879
880	c-----------------------------------------------------------------------
881	c Variable 2d :
882	c-----------------------------------------------------------------------
883
884	c Relief
885	! topoflag = F: we keep the existing topography
886	! = T: we read it from the Relief.nc file
887	! topoflag need to be in the specific run.def for newstart
888	topoflag = . FALSE .
889	CALL getin('topoflag',topoflag)
890	! notopo = T: we go back to flat surface
891	notopo = .FALSE.
892	CALL getin('notopo',notopo)
893
894	print*,zmeaold(2,1:10)
895
896	IF ( topoflag ) then
897	print*,"Topography (phis) read in file Relief.nc"
898	print*,"For Venus, map directly inverted in Relief.nc"
899	CALL startget_phys2d('surfgeo',iip1,jjp1,rlonv,rlatu,phis,0.0,
900	. jjm ,rlonu,rlatv,.true.)
901	c needed ?
902	phis(iip1,:) = phis(1,:)
903
904	CALL startget_phys1d('zmea',iip1,jjp1,rlonv,rlatu,ngridmx,zmea,
905	. 0.0,jjm,rlonu,rlatv,.true.)
906	CALL startget_phys1d('zstd',iip1,jjp1,rlonv,rlatu,ngridmx,zstd,
907	. 0.0,jjm,rlonu,rlatv,.true.)
908	CALL startget_phys1d('zsig',iip1,jjp1,rlonv,rlatu,ngridmx,zsig,
909	. 0.0,jjm,rlonu,rlatv,.true.)
910	CALL startget_phys1d('zgam',iip1,jjp1,rlonv,rlatu,ngridmx,zgam,
911	. 0.0,jjm,rlonu,rlatv,.true.)
912	CALL startget_phys1d('zthe',iip1,jjp1,rlonv,rlatu,ngridmx,zthe,
913	. 0.0,jjm,rlonu,rlatv,.true.)
914	CALL startget_phys1d('zpic',iip1,jjp1,rlonv,rlatu,ngridmx,zpic,
915	. 0.0,jjm,rlonu,rlatv,.true.)
916	CALL startget_phys1d('zval',iip1,jjp1,rlonv,rlatu,ngridmx,zval,
917	. 0.0,jjm,rlonu,rlatv,.true.)
918
919	ELSE IF ( notopo ) THEN
920	print*,'Flattening the topography'
921	phis=0.
922	zmea=0.
923	zstd=0.
924	zsig=0.
925	zgam=0.
926	zthe=0.
927	zpic=0.
928	zval=0.
929	ELSE
930	print*,'Using existing topography'
931	call interp_horiz (phisold,phis,imold,jmold,iim,jjm,1,
932	& rlonuold,rlatvold,rlonu,rlatv)
933
934	call interp_horiz (zmeaold,zmeaS,imold,jmold,iim,jjm,1,
935	& rlonuold,rlatvold,rlonu,rlatv)
936	call gr_dyn_fi (1,iip1,jjp1,ngridmx,zmeaS,zmea)
937	call interp_horiz (zstdold,zstdS,imold,jmold,iim,jjm,1,
938	& rlonuold,rlatvold,rlonu,rlatv)
939	call gr_dyn_fi (1,iip1,jjp1,ngridmx,zstdS,zstd)
940	call interp_horiz (zsigold,zsigS,imold,jmold,iim,jjm,1,
941	& rlonuold,rlatvold,rlonu,rlatv)
942	call gr_dyn_fi (1,iip1,jjp1,ngridmx,zsigS,zsig)
943	call interp_horiz (zgamold,zgamS,imold,jmold,iim,jjm,1,
944	& rlonuold,rlatvold,rlonu,rlatv)
945	call gr_dyn_fi (1,iip1,jjp1,ngridmx,zgamS,zgam)
946	call interp_horiz (ztheold,ztheS,imold,jmold,iim,jjm,1,
947	& rlonuold,rlatvold,rlonu,rlatv)
948	call gr_dyn_fi (1,iip1,jjp1,ngridmx,ztheS,zthe)
949	call interp_horiz (zpicold,zpicS,imold,jmold,iim,jjm,1,
950	& rlonuold,rlatvold,rlonu,rlatv)
951	call gr_dyn_fi (1,iip1,jjp1,ngridmx,zpicS,zpic)
952	call interp_horiz (zvalold,zvalS,imold,jmold,iim,jjm,1,
953	& rlonuold,rlatvold,rlonu,rlatv)
954	call gr_dyn_fi (1,iip1,jjp1,ngridmx,zvalS,zval)
955	ENDIF
956
957	print*,"New surface geopotential OK"
958
959	c Lat et lon pour physique
960	do i=1,iip1
961	rlatS(i,:)=rlatu(:)*180./pi
962	enddo
963	call gr_dyn_fi (1,iip1,jjp1,ngridmx,rlatS,rlat)
964
965	do j=2,jjm
966	rlonS(:,j)=rlonv(:)*180./pi
967	enddo
968	rlonS(:,1)=0.
969	rlonS(:,jjp1)=0.
970	call gr_dyn_fi (1,iip1,jjp1,ngridmx,rlonS,rlon)
971
972	c Temperature de surface
973	! razTS need to be in the specific run.def for newstart
974	razTS = . FALSE .
975	CALL getin('razTS',razTS)
976
977	if (razTS) then
978	tsurf(:) = 735.
979	else
980	call interp_horiz (tsurfold,tsurfS,imold,jmold,iim,jjm,1,
981	& rlonuold,rlatvold,rlonu,rlatv)
982	call gr_dyn_fi (1,iip1,jjp1,ngridmx,tsurfS,tsurf)
983	c write(44,*) 'tsurf', tsurf
984	endif
985
986	c Temperature du sous-sol
987	if (razTS) then
988	tsoil(:,:)=735.
989	else
990	call interp_horiz(tsoilold,tsoilS,
991	& imold,jmold,iim,jjm,nsoilmx,
992	& rlonuold,rlatvold,rlonu,rlatv)
993	call gr_dyn_fi (nsoilmx,iip1,jjp1,ngridmx,tsoilS,tsoil)
994	c write(45,*) 'tsoil',tsoil
995	endif
996
997	! CHANGING ALBEDO: may be done through run.def
998	albedoflag = . FALSE .
999	CALL getin('albedoflag',albedoflag)
1000
1001	IF ( albedoflag ) then
1002	print*,"Albedo is reinitialized to the albedo value in run.def"
1003	print*,"We may want to consider a map later on..."
1004	albedo=0.1
1005	CALL getin('albedo',albedo)
1006	albe=albedo
1007	ELSE
1008	call interp_horiz (albeold,albeS,imold,jmold,iim,jjm,1,
1009	& rlonuold,rlatvold,rlonu,rlatv)
1010	call gr_dyn_fi (1,iip1,jjp1,ngridmx,albeS,albe)
1011	ENDIF
1012
1013	call interp_horiz (radsolold,radsolS,imold,jmold,iim,jjm,1,
1014	& rlonuold,rlatvold,rlonu,rlatv)
1015	call gr_dyn_fi (1,iip1,jjp1,ngridmx,radsolS,radsol)
1016
1017	call interp_horiz (sollwold,sollwS,imold,jmold,iim,jjm,1,
1018	& rlonuold,rlatvold,rlonu,rlatv)
1019	call gr_dyn_fi (1,iip1,jjp1,ngridmx,sollwS,sollw)
1020
1021	call interp_horiz (solswold,solswS,imold,jmold,iim,jjm,1,
1022	& rlonuold,rlatvold,rlonu,rlatv)
1023	call gr_dyn_fi (1,iip1,jjp1,ngridmx,solswS,solsw)
1024
1025	call interp_horiz (fderold,fderS,imold,jmold,iim,jjm,1,
1026	& rlonuold,rlatvold,rlonu,rlatv)
1027	call gr_dyn_fi (1,iip1,jjp1,ngridmx,fderS,fder)
1028
1029	call interp_horiz (dlwold,dlwS,imold,jmold,iim,jjm,1,
1030	& rlonuold,rlatvold,rlonu,rlatv)
1031	call gr_dyn_fi (1,iip1,jjp1,ngridmx,dlwS,dlw)
1032
1033	call interp_horiz (sollwdownold,sollwdownS,imold,jmold,iim,jjm,1,
1034	& rlonuold,rlatvold,rlonu,rlatv)
1035	call gr_dyn_fi (1,iip1,jjp1,ngridmx,sollwdownS,sollwdown)
1036
1037	print*,"Nouvelles var physiques OK"
1038
1039	c-----------------------------------------------------------------------
1040	c Traitement special de la pression au sol :
1041	c-----------------------------------------------------------------------
1042
1043	c Extrapolation la pression dans la nouvelle grille
1044	call interp_horiz(psold,ps,imold,jmold,iim,jjm,1,
1045	& rlonuold,rlatvold,rlonu,rlatv)
1046
1047	c On assure la conservation de la masse de l'atmosphere
1048	c--------------------------------------------------------------
1049
1050	ptotal = 0.
1051	DO j=1,jjp1
1052	DO i=1,iim
1053	ptotal=ptotal+ps(i,j)*aire(i,j)/g
1054	END DO
1055	END DO
1056
1057	write(,)
1058	write(,)'Ancienne grille: masse de l atm :',ptotalold
1059	write(,)'Nouvelle grille: masse de l atm :',ptotal
1060	write (,) 'Ratio new atm./ old atm =', ptotal/ptotalold
1061	write(,)
1062
1063
1064	DO j=1,jjp1
1065	DO i=1,iip1
1066	ps(i,j)=ps(i,j) * ptotalold/ptotal
1067	END DO
1068	END DO
1069
1070	c la pression de surface et les temperatures ne sont pas reequilibrees en fonction
1071	c de la nouvelle topographie...
1072	c Si l'ajustement inevitable du debut pose des problemes, voir le newstart martien.
1073
1074	print*,"Nouvelle ps OK"
1075	print*,"If changes done on topography, beware !"
1076	print*,"Some time may be needed for adjustments at the beginning"
1077	print*,"so if unstable, relax the timestep and/or dissipation."
1078
1079	c-----------------------------------------------------------------------
1080	c Variable 3d :
1081	c-----------------------------------------------------------------------
1082
1083	CALL pression(ip1jmp1, ap, bp, ps, p3d)
1084	if (disvert_type==1) then
1085	CALL exner_hyb( ip1jmp1, ps, p3d,pks, pk, pkf )
1086	else ! we assume that we are in the disvert_type==2 case
1087	CALL exner_milieu( ip1jmp1, ps, p3d, pks, pk, pkf )
1088	endif
1089
1090	c temperatures atmospheriques
1091
1092	c ATTENTION: peut servir, mais bon...
1093	c modif: profil uniforme
1094	c do l=1,lmold
1095	c do j=1,jmold+1
1096	c do i=1,imold+1
1097	c told(i,j,l)=told(1,jmold/2,l)
1098	c enddo
1099	c enddo
1100	c enddo
1101
1102	! raztemp need to be in the specific run.def for newstart
1103	raztemp = . FALSE .
1104	CALL getin('raztemp',raztemp)
1105
1106	! Reinitialisation of temperature to VIRA profile lisse
1107	if (raztemp) then
1108
1109	allocate(tvira(0:lmold),dzst(0:lmold),zkm(0:lmold))
1110	print*,"Venus = temperature initiale imposee = VIRA lisse "
1111	dzst(0) = 0.0
1112	dzst(1) = -log(p3d(1,1,2)/preff)*r/g
1113	do l=2,lmold
1114	dzst(l)=-log(p3d(1,1,l+1)/p3d(1,1,l))*r/g
1115	enddo
1116	tvira(0) = 735.
1117	zkm(0) = 0.0
1118	do l=1,lmold
1119	zkm(l) = zkm(l-1)+tvira(l-1)*dzst(l)/1000. ! approx avec T(l-1)
1120	if(zkm(l).lt.60.) then
1121	tvira(l)=735.-7.95*zkm(l)
1122	else
1123	tvira(l)=AMAX1(258.-3.*(zkm(l)-60.),168.)
1124	endif
1125	zkm(l) = zkm(l-1)+(tvira(l-1)+tvira(l))/2.*dzst(l)/1000.
1126	enddo
1127	do l=1,lmold
1128	do j=1,jmold+1
1129	do i=1,imold+1
1130	told(i,j,l)=tvira(l)
1131	enddo
1132	enddo
1133	enddo
1134	endif ! end raztemp
1135
1136	write (,) 'told ', told (1,jmold+1,1) ! INFO
1137	call interp_vert
1138	& (told,var,lmold,llm,apold,bpold,ap,bp,
1139	& psold,(imold+1)*(jmold+1))
1140	write (,) 'var ', var (1,jmold+1,1) ! INFO
1141	call interp_horiz(var,t,imold,jmold,iim,jjm,llm,
1142	& rlonuold,rlatvold,rlonu,rlatv)
1143	write (,) 'T ', T(1,jjp1,1) ! INFO
1144	! pour info:
1145	! Si extension verticale, la T est extrapolee constante au-dessus de lmold
1146
1147	c passage grille physique pour restartphy.nc
1148	call gr_dyn_fi (llm,iip1,jjp1,ngridmx,T,t_fi)
1149
1150	! ADAPTATION GCM POUR CP(T)
1151	c passage en temperature potentielle
1152	call t2tpot(ip1jmp1*llm,T,teta,pk)
1153	c on assure la periodicite
1154	teta(iip1,:,:) = teta(1,:,:)
1155
1156	! RESETING U TO uini: may be done through run.def
1157	uini = .FALSE.
1158	CALL getin('uini',uini)
1159	! RESETING U TO 0: may be done through run.def
1160	razvitu = . FALSE .
1161	CALL getin('razvitu',razvitu)
1162	razvitv = . FALSE .
1163	CALL getin('razvitv',razvitv)
1164
1165	c calcul des champ de vent; passage en vent covariant
1166	write (,) 'uold ', uold (1,2,1) ! INFO
1167	call interp_vert
1168	& (uold,var,lmold,llm,apold,bpold,ap,bp,
1169	& psold,(imold+1)*(jmold+1))
1170	write (,) 'var ', var (1,2,1) ! INFO
1171	call interp_horiz(var,us,imold,jmold,iim,jjm,llm,
1172	& rlonuold,rlatvold,rlonu,rlatv)
1173	write (,) 'us ', us (1,2,1) ! INFO
1174
1175	call interp_vert
1176	& (vold,var,lmold,llm,
1177	& apold,bpold,ap,bp,psold,(imold+1)*(jmold+1))
1178	call interp_horiz(var,vs,imold,jmold,iim,jjm,llm,
1179	& rlonuold,rlatvold,rlonu,rlatv)
1180	call scal_wind(us,vs,unat,vnat)
1181	! Reseting u=0
1182	if ((razvitu).and..not.(uini)) then
1183	unat(:,:,:) = 0.
1184	endif
1185	! Reseting u=uini
1186	if ((uini).and..not.(razvitu)) then
1187	do j=1,jjp1
1188	do l=1,llm
1189	if (p3d(1,j,l).gt.3e3) then
1190	unat(:,j,l) = -110./8.03*log(p3d(:,j,l)/9.2e6)
1191	else
1192	unat(:,j,l) = 110./6.62*(log(p3d(:,j,l)/9.2e6)+14.65)
1193	endif
1194	if (abs(rlatS(1,j)).gt.50.) then
1195	unat(:,j,l)=unat(:,j,l)*(90.-abs(rlatS(:,j)))/40.
1196	endif
1197	enddo
1198	enddo
1199	endif
1200	! incompatible options
1201	if ((uini).and.(razvitu)) then
1202	print*,"You have to choose between razvitu and uini..."
1203	stop
1204	endif
1205	write (,) 'unat ', unat (1,2,1) ! INFO
1206
1207	do l=1,llm
1208	do j = 1, jjp1
1209	do i=1,iip1
1210	ucov( i,j,l ) = unat( i,j,l ) * cu(i,j)
1211	! pour info:
1212	! Si extension verticale, on impose u=0 au-dessus de lmold
1213	if (l.gt.lmold) ucov( i,j,l ) = 0
1214	end do
1215	end do
1216	end do
1217	write (,) 'ucov ', ucov (1,2,1) ! INFO
1218	c write(48,*) 'ucov',ucov
1219	! Reseting v=0
1220	if (razvitv) then
1221	vnat(:,:,:) = 0.
1222	endif
1223	write (,) 'vnat ', vnat (1,2,1) ! INFO
1224	do l=1,llm
1225	do j = 1, jjm
1226	do i=1,iim
1227	vcov( i,j,l ) = vnat( i,j,l ) * cv(i,j)
1228	! pour info:
1229	! Si extension verticale, on impose v=0 au-dessus de lmold
1230	if (l.gt.lmold) vcov( i,j,l ) = 0
1231	end do
1232	vcov( iip1,j,l ) = vcov( 1,j,l )
1233	end do
1234	end do
1235	c write(49,*) 'ucov',vcov
1236
1237	c masse: on la recalcule (ps a été ajustée pour conserver la masse totale)
1238	call massdair(p3d,masse)
1239
1240	c traceurs 3D
1241	do iq = 1, nqtot
1242	call interp_vert(qold(1,1,1,iq),var,lmold,llm,
1243	& apold,bpold,ap,bp,psold,(imold+1)*(jmold+1))
1244	call interp_horiz(var,q(1,1,1,iq),imold,jmold,iim,jjm,llm,
1245	& rlonuold,rlatvold,rlonu,rlatv)
1246	enddo
1247
1248	print*,"Nouvelles var dynamiques OK"
1249
1250	c=======================================================================
1251	c Ecriture des restart.nc et restartphy.nc :
1252	c=======================================================================
1253
1254	call writerestart('restart.nc',tab_cntrl_dyn,
1255	. phis,vcov,ucov,teta,masse,q,ps)
1256
1257	print*,"restart OK"
1258
1259	call writerestartphy('restartphy.nc',tab_cntrl_fi,ngridmx,llm,
1260	. rlat,rlon,tsurf,tsoil,albe,solsw, sollw,fder,dlw,
1261	. sollwdown,radsol,
1262	. zmea,zstd,zsig,zgam,zthe,zpic,zval,
1263	. t_fi)
1264
1265	print*,"restartphy OK"
1266
1267	end

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

Download in other formats:

Original Format