Context Navigation

newstart.F @ 1243

Last change on this file since 1243 was 1056, checked in by slebonnois, 11 years ago
SL: Titan runs ! see DOC/chantiers/commit_importants.log
File size: 37.1 KB

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

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

Download in other formats:

Original Format