Context Navigation

lect_start_archive.F @ 1036

Last change on this file since 1036 was 1036, checked in by emillour, 12 years ago

Mars GCM: (a first step towards using parallel dynamics)

IMPORTANT CHANGE: Implemented dynamic tracers. It is no longer necessary to compile the model with the '-t #' option, number of tracers is simply read from tracer.def file (as before). Adapted makegcm_* scripts (and co.) accordingly. Technical aspects of the switch to dynamic tracers are:
- advtrac.h (in dyn3d) removed and replaced by module infotrac.F
- tracer.h (in phymars) removed and replaced by module tracer_mod.F90 (which contains nqmx, the number of tracers, etc. and can be used anywhere in the physics).
Included some side cleanups: removed unused files (in dyn3d) anldoppler2.F, anl_mcdstats.F and anl_stats-diag.F, and all the unecessary dimensions.* files in grid/dimension.
Checked that changes are clean and that GCM yields identical results (in debug mode) to previous svn version.

File size: 45.0 KB

Rev	Line
[1036]	1	SUBROUTINE lect_start_archive(nqtot,date,tsurf,tsoil,emis,q2,
[38]	2	& t,ucov,vcov,ps,co2ice,h,phisold_newgrid,
	3	& q,qsurf,surfith,nid)
	4	c=======================================================================
	5	c
	6	c
	7	c Auteur: 05/1997 , 12/2003 : coord hybride FF
	8	c ------
	9	c
	10	c
	11	c Objet: Lecture des variables d'un fichier "start_archive"
	12	c Plus besoin de régler ancienne valeurs grace
	13	c a l'allocation dynamique de memoire (Yann Wanherdrick)
	14	c
	15	c
	16	c
	17	c=======================================================================
[1036]	18	use infotrac, only: tnom
[38]	19	implicit none
	20
	21	#include "dimensions.h"
	22	#include "dimphys.h"
	23	#include "surfdat.h"
	24	#include "comsoil.h"
	25	#include "dimradmars.h"
	26	#include "yomaer.h"
	27	#include "planete.h"
	28	#include "paramet.h"
	29	#include "comconst.h"
	30	#include "comvert.h"
	31	#include "comgeom2.h"
	32	#include "control.h"
	33	#include "logic.h"
	34	#include "description.h"
	35	#include "ener.h"
	36	#include "temps.h"
	37	#include "lmdstd.h"
	38	#include "netcdf.inc"
	39	!#include "tracer.h"
[1036]	40	!#include"advtrac.h"
[38]	41	c=======================================================================
	42	c Declarations
	43	c=======================================================================
	44
	45	c Old variables dimensions (from file)
	46	c------------------------------------
	47	INTEGER imold,jmold,lmold,nsoilold,nqold
	48
	49	c et autres:
	50	c----------
[1036]	51	integer,intent(in) :: nqtot ! number of advected tracers
[38]	52
	53	c Variables pour les lectures des fichiers "ini"
	54	c--------------------------------------------------
	55	! INTEGER sizei,
	56	integer timelen,dimid
	57	! INTEGER length
	58	! parameter (length = 100)
	59	INTEGER tab0
	60	INTEGER isoil,iq,iqmax
	61	CHARACTER*2 str2
	62
	63	! REAL dimfirst(4) ! tableau contenant les 1ers elements des dimensions
	64
	65	! REAL dimlast(4) ! tableau contenant les derniers elements des dimensions
	66
	67	! REAL dimcycl(4) ! tableau contenant les periodes des dimensions
	68	! CHARACTER*120 dimsource
	69	! CHARACTER*16 dimname
	70	! CHARACTER*80 dimtitle
	71	! CHARACTER*40 dimunits
	72	! INTEGER dimtype
	73
	74	! INTEGER dimord(4) ! tableau contenant l''ordre
	75	! data dimord /1,2,3,4/ ! de sortie des dimensions
	76
	77	! INTEGER vardim(4)
	78	REAL date
	79	INTEGER memo
	80	! character (len=50) :: tmpname
	81
	82	c Variable histoire
	83	c------------------
	84	REAL vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm) ! vents covariants
	85	REAL h(iip1,jjp1,llm),ps(iip1,jjp1)
[1036]	86	REAL q(iip1,jjp1,llm,nqtot),qtot(iip1,jjp1,llm)
[38]	87
	88	c autre variables dynamique nouvelle grille
	89	c------------------------------------------
	90
	91	c!-*-
	92	! integer klatdat,klongdat
	93	! PARAMETER (klatdat=180,klongdat=360)
	94
	95	c Physique sur grille scalaire
	96	c----------------------------
	97
	98	c variable physique
	99	c------------------
	100	REAL tsurf(ngridmx) ! surface temperature
	101	REAL tsoil(ngridmx,nsoilmx) ! soil temperature
	102	REAL co2ice(ngridmx) ! CO2 ice layer
	103	REAL emis(ngridmx)
[1036]	104	REAL q2(ngridmx,nlayermx+1),qsurf(ngridmx,nqtot)
[38]	105	c REAL phisfi(ngridmx)
	106
	107	INTEGER i,j,l
	108	INTEGER nid,nvarid
	109	c REAL year_day,periheli,aphelie,peri_day
	110	c REAL obliquit,z0,emin_turb,lmixmin
	111	c REAL emissiv,emisice(2),albedice(2),tauvis
	112	c REAL iceradius(2) , dtemisice(2)
	113
	114	! EXTERNAL RAN1
	115	! REAL RAN1
	116	! EXTERNAL geopot,inigeom
	117	integer ierr
	118	! integer ismin
	119	! external ismin
	120	! CHARACTER*80 datapath
	121	integer, dimension(4) :: start,count
	122
	123	c Variable nouvelle grille naturelle au point scalaire
	124	c------------------------------------------------------
	125	real us(iip1,jjp1,llm),vs(iip1,jjp1,llm)
	126	REAL phisold_newgrid(iip1,jjp1)
	127	REAL t(iip1,jjp1,llm)
	128	real tsurfS(iip1,jjp1),tsoilS(iip1,jjp1,nsoilmx)
	129	real inertiedatS(iip1,jjp1,nsoilmx)
	130	real co2iceS(iip1,jjp1),emisS(iip1,jjp1)
[1036]	131	REAL q2S(iip1,jjp1,llm+1),qsurfS(iip1,jjp1,nqtot)
[38]	132
	133	real ptotal, co2icetotal
	134
	135	c Var intermediaires : vent naturel, mais pas coord scalaire
	136	c-----------------------------------------------------------
	137	real vnat(iip1,jjm,llm),unat(iip1,jjp1,llm)
	138
	139
	140	c Variable de l'ancienne grille
	141	c---------------------------------------------------------
	142
	143	real, dimension(:), allocatable :: timelist
	144	real, dimension(:), allocatable :: rlonuold, rlatvold
	145	real, dimension(:), allocatable :: rlonvold, rlatuold
	146	real, dimension(:), allocatable :: apsold,bpsold
	147	real, dimension(:), allocatable :: mlayerold
	148	real, dimension(:,:,:), allocatable :: uold,vold,told,q2old
	149	real, dimension(:,:,:), allocatable :: tsoilold,qsurfold
	150	real, dimension(:,:,:),allocatable :: tsoiloldnew
	151	! tsoiloldnew: old soil values, but along new subterranean grid
	152	real, dimension(:,:,:), allocatable :: inertiedatold
	153	! inertiedatoldnew: old inertia values, but along new subterranean grid
	154	real, dimension(:,:,:), allocatable :: inertiedatoldnew
	155	real, dimension(:,:), allocatable :: psold,phisold
	156	real, dimension(:,:), allocatable :: co2iceold,tsurfold
	157	real, dimension(:,:), allocatable :: emisold
	158	real, dimension(:,:,:,:), allocatable :: qold
	159
	160	real tab_cntrl(100)
	161
	162	real ptotalold, co2icetotalold
	163
	164	logical :: olddepthdef=.false. ! flag
	165	! olddepthdef=.true. if soil depths are in 'old' (unspecified) format
	166	logical :: depthinterpol=.false. ! flag
	167	! depthinterpol=.true. if interpolation will be requiered
	168	logical :: therminertia_3D=.true. ! flag
	169	! therminertia_3D=.true. if thermal inertia is 3D and read from datafile
	170	c Variable intermediaires iutilise pour l'extrapolation verticale
	171	c----------------------------------------------------------------
	172	real, dimension(:,:,:), allocatable :: var,varp1
	173	real, dimension(:), allocatable :: oldgrid, oldval
	174	real, dimension(:), allocatable :: newval
	175	! real, dimension(:), allocatable :: oldmlayer
	176
	177	real surfith(iip1,jjp1) ! surface thermal inertia
	178	! real surfithfi(ngridmx)
	179	! surface thermal inertia at old horizontal grid resolution
	180	real, dimension(:,:), allocatable :: surfithold
	181
	182	! flag which identifies if archive file is using old tracer names (qsurf01,...)
	183	logical :: oldtracernames=.false.
	184	integer :: counter
	185	character(len=30) :: txt ! to store some text
[146]	186	real :: tmpval ! to store a temporary variable/value
[38]	187
	188	c=======================================================================
	189
	190	! 0. Preliminary stuff
	191
	192	! check if tracers follow old naming convention (q01, q02, q03, ...)
	193	counter=0
[1036]	194	do iq=1,nqtot
[38]	195	txt= " "
	196	write(txt,'(a1,i2.2)')'q',iq
	197	ierr=NF_INQ_VARID(nid,txt,nvarid)
	198	if (ierr.ne.NF_NOERR) then
	199	! did not find old tracer name
	200	exit ! might as well stop here
	201	else
	202	! found old tracer name
	203	counter=counter+1
	204	endif
	205	enddo
[1036]	206	if (counter.eq.nqtot) then
[38]	207	write(,) "lect_start_archive: tracers seem to follow old ",
	208	& "naming convention (q01, q02,...)"
	209	oldtracernames=.true.
	210	endif
	211
	212
	213	!-----------------------------------------------------------------------
	214	! 1. Read data dimensions (i.e. size and length)
	215	!-----------------------------------------------------------------------
	216
	217	! 1.2 Read the various dimension lengths of data in file
	218
	219	ierr= NF_INQ_DIMID(nid,"Time",dimid)
	220	if (ierr.ne.NF_NOERR) then
	221	ierr= NF_INQ_DIMID(nid,"temps",dimid)
	222	endif
	223	ierr= NF_INQ_DIMLEN(nid,dimid,timelen)
	224
	225	ierr= NF_INQ_DIMID(nid,"latitude",dimid)
	226	if (ierr.ne.NF_NOERR) then
	227	ierr= NF_INQ_DIMID(nid,"rlatu",dimid)
	228	endif
	229	ierr= NF_INQ_DIMLEN(nid,dimid,jmold)
	230	jmold=jmold-1
	231
	232	ierr= NF_INQ_DIMID(nid,"longitude",dimid)
	233	if (ierr.ne.NF_NOERR) then
	234	ierr= NF_INQ_DIMID(nid,"rlonv",dimid)
	235	endif
	236	ierr= NF_INQ_DIMLEN(nid,dimid,imold)
	237	imold=imold-1
	238
	239	ierr= NF_INQ_DIMID(nid,"altitude",dimid)
	240	if (ierr.ne.NF_NOERR) then
	241	ierr= NF_INQ_DIMID(nid,"sig_s",dimid)
	242	endif
	243	ierr= NF_INQ_DIMLEN(nid,dimid,lmold)
	244
	245	nqold=0
	246	do
	247	write(str2,'(i2.2)') nqold+1
	248	ierr= NF_INQ_VARID(nid,'q'//str2,dimid)
	249	! write(,) 'q'//str2
	250	if (ierr.eq.NF_NOERR) then
	251	nqold=nqold+1
	252	else
	253	exit
	254	endif
	255	enddo
	256
	257	! 1.2.1 find out the # of subsurface_layers
	258	nsoilold=0 !dummy initialisation
	259	ierr=NF_INQ_DIMID(nid,"subsurface_layers",dimid)
	260	if (ierr.eq.NF_NOERR) then
	261	ierr=NF_INQ_DIMLEN(nid,dimid,nsoilold)
	262	if (ierr.ne.NF_NOERR) then
	263	write(,)'lec_start_archive: ',
	264	& 'Failed reading subsurface_layers length'
	265	endif
	266	else
	267	write(,)"lec_start_archive: did not find subsurface_layers"
	268	endif
	269
	270	if (nsoilold.eq.0) then ! 'old' archive format;
	271	! must use Tg//str2 fields to compute nsoilold
	272	write(,)"lec_start_archive: building nsoilold from Tg fields"
	273	do
	274	write(str2,'(i2.2)') nsoilold+1
	275	ierr=NF_INQ_VARID(nid,'Tg'//str2,dimid)
	276	if (ierr.eq.NF_NOERR) then
	277	nsoilold=nsoilold+1
	278	else
	279	exit
	280	endif
	281	enddo
	282	endif
	283
	284
	285	if (nsoilold.ne.nsoilmx) then ! interpolation will be required
	286	depthinterpol=.true.
	287	endif
	288
	289	! 1.3 Report dimensions
	290
	291	write(,) "Start_archive dimensions:"
	292	write(,) "longitude: ",imold
	293	write(,) "latitude: ",jmold
	294	write(,) "altitude: ",lmold
	295	write(,) "tracers: ",nqold
	296	write(,) "subsurface_layers: ",nsoilold
	297	if (depthinterpol) then
	298	write(,) " => Warning, nsoilmx= ",nsoilmx
	299	write(,) ' which implies that you want subterranean interpola
	300	&tion.'
	301	write(,) ' Otherwise, set nsoilmx -in dimphys.h- to: ',nsoilold
	302	endif
	303	write(,) "time lenght: ",timelen
	304	write(,)
	305
	306	!-----------------------------------------------------------------------
	307	! 2. Allocate arrays to store datasets
	308	!-----------------------------------------------------------------------
	309
	310	allocate(timelist(timelen))
	311	allocate(rlonuold(imold+1), rlatvold(jmold))
	312	allocate(rlonvold(imold+1), rlatuold(jmold+1))
	313	allocate (apsold(lmold),bpsold(lmold))
	314	allocate(uold(imold+1,jmold+1,lmold))
	315	allocate(vold(imold+1,jmold+1,lmold))
	316	allocate(told(imold+1,jmold+1,lmold))
	317	allocate(psold(imold+1,jmold+1))
	318	allocate(phisold(imold+1,jmold+1))
[1036]	319	allocate(qold(imold+1,jmold+1,lmold,nqtot))
[38]	320	allocate(co2iceold(imold+1,jmold+1))
	321	allocate(tsurfold(imold+1,jmold+1))
	322	allocate(emisold(imold+1,jmold+1))
	323	allocate(q2old(imold+1,jmold+1,lmold+1))
	324	! allocate(tsoilold(imold+1,jmold+1,nsoilmx))
	325	allocate(tsoilold(imold+1,jmold+1,nsoilold))
	326	allocate(tsoiloldnew(imold+1,jmold+1,nsoilmx))
	327	allocate(inertiedatold(imold+1,jmold+1,nsoilold)) ! soil thermal inertia
	328	allocate(inertiedatoldnew(imold+1,jmold+1,nsoilmx))
	329	! surface thermal inertia at old horizontal grid resolution
	330	allocate(surfithold(imold+1,jmold+1))
	331	allocate(mlayerold(nsoilold))
[1036]	332	allocate(qsurfold(imold+1,jmold+1,nqtot))
[38]	333
	334	allocate(var (imold+1,jmold+1,llm))
	335	allocate(varp1 (imold+1,jmold+1,llm+1))
	336
	337	write(,) 'q2',ngridmx,nlayermx+1
	338	write(,) 'q2S',iip1,jjp1,llm+1
	339	write(,) 'q2old',imold+1,jmold+1,lmold+1
	340
	341	!-----------------------------------------------------------------------
	342	! 3. Read time-independent data
	343	!-----------------------------------------------------------------------
	344
	345	C-----------------------------------------------------------------------
	346	c 3.1. Lecture du tableau des parametres du run
	347	c (pour la lecture ulterieure de "ptotalold" et "co2icetotalold")
	348	c-----------------------------------------------------------------------
	349	c
	350	ierr = NF_INQ_VARID (nid, "controle", nvarid)
	351	IF (ierr .NE. NF_NOERR) THEN
[146]	352	PRINT*, "Lect_start_archive: <controle> is missing"
[38]	353	CALL abort
	354	ENDIF
	355	#ifdef NC_DOUBLE
	356	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tab_cntrl)
	357	#else
	358	ierr = NF_GET_VAR_REAL(nid, nvarid, tab_cntrl)
	359	#endif
	360	IF (ierr .NE. NF_NOERR) THEN
[146]	361	PRINT*, "lect_start_archive: Failed loading <controle>"
[38]	362	CALL abort
	363	ENDIF
	364	c
	365	tab0 = 50
	366
	367	c-----------------------------------------------------------------------
	368	c 3.2 Lecture des longitudes et latitudes
	369	c-----------------------------------------------------------------------
	370	c
	371	ierr = NF_INQ_VARID (nid, "rlonv", nvarid)
	372	IF (ierr .NE. NF_NOERR) THEN
[146]	373	PRINT*, "lect_start_archive: <rlonv> is missing"
[38]	374	CALL abort
	375	ENDIF
	376	#ifdef NC_DOUBLE
	377	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlonvold)
	378	#else
	379	ierr = NF_GET_VAR_REAL(nid, nvarid, rlonvold)
	380	#endif
	381	IF (ierr .NE. NF_NOERR) THEN
[146]	382	PRINT*, "lect_start_archive: Failed loading <rlonv>"
[38]	383	CALL abort
	384	ENDIF
	385	c
	386	ierr = NF_INQ_VARID (nid, "rlatu", nvarid)
	387	IF (ierr .NE. NF_NOERR) THEN
[146]	388	PRINT*, "lect_start_archive: <rlatu> is missing"
[38]	389	CALL abort
	390	ENDIF
	391	#ifdef NC_DOUBLE
	392	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlatuold)
	393	#else
	394	ierr = NF_GET_VAR_REAL(nid, nvarid, rlatuold)
	395	#endif
	396	IF (ierr .NE. NF_NOERR) THEN
[146]	397	PRINT*, "lect_start_archive: Failed loading <rlatu>"
[38]	398	CALL abort
	399	ENDIF
	400	c
	401	ierr = NF_INQ_VARID (nid, "rlonu", nvarid)
	402	IF (ierr .NE. NF_NOERR) THEN
[146]	403	PRINT*, "lect_start_archive: <rlonu> is missing"
[38]	404	CALL abort
	405	ENDIF
	406	#ifdef NC_DOUBLE
	407	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlonuold)
	408	#else
	409	ierr = NF_GET_VAR_REAL(nid, nvarid, rlonuold)
	410	#endif
	411	IF (ierr .NE. NF_NOERR) THEN
[146]	412	PRINT*, "lect_start_archive: Failed loading <rlonu>"
[38]	413	CALL abort
	414	ENDIF
	415	c
	416	ierr = NF_INQ_VARID (nid, "rlatv", nvarid)
	417	IF (ierr .NE. NF_NOERR) THEN
[146]	418	PRINT*, "lect_start_archive: <rlatv> is missing"
[38]	419	CALL abort
	420	ENDIF
	421	#ifdef NC_DOUBLE
	422	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlatvold)
	423	#else
	424	ierr = NF_GET_VAR_REAL(nid, nvarid, rlatvold)
	425	#endif
	426	IF (ierr .NE. NF_NOERR) THEN
[146]	427	PRINT*, "lect_start_archive: Failed loading <rlatv>"
[38]	428	CALL abort
	429	ENDIF
	430	c
	431
	432	c-----------------------------------------------------------------------
	433	c 3.3. Lecture des niveaux verticaux
	434	c-----------------------------------------------------------------------
	435	c
	436	ierr = NF_INQ_VARID (nid, "aps", nvarid)
	437	IF (ierr .NE. NF_NOERR) THEN
[146]	438	PRINT*, "lect_start_archive: <aps> is missing"
[38]	439	apsold=0
	440	PRINT*, "<aps> set to 0"
	441	ELSE
	442	#ifdef NC_DOUBLE
	443	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, apsold)
	444	#else
	445	ierr = NF_GET_VAR_REAL(nid, nvarid, apsold)
	446	#endif
	447	IF (ierr .NE. NF_NOERR) THEN
[146]	448	PRINT*, "lect_start_archive: Failed loading <aps>"
[38]	449	ENDIF
	450	ENDIF
	451	c
	452	ierr = NF_INQ_VARID (nid, "bps", nvarid)
	453	IF (ierr .NE. NF_NOERR) THEN
[146]	454	PRINT*, "lect_start_archive: <bps> is missing"
[38]	455	PRINT*, "It must be an old start_archive, lets look for sig_s"
	456	ierr = NF_INQ_VARID (nid, "sig_s", nvarid)
	457	IF (ierr .NE. NF_NOERR) THEN
	458	PRINT*, "Nothing to do..."
	459	CALL abort
	460	ENDIF
	461	ENDIF
	462	#ifdef NC_DOUBLE
	463	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, bpsold)
	464	#else
	465	ierr = NF_GET_VAR_REAL(nid, nvarid, bpsold)
	466	#endif
	467	IF (ierr .NE. NF_NOERR) THEN
[146]	468	PRINT*, "lect_start_archive: Failed loading <bps>"
[38]	469	CALL abort
	470	END IF
	471
	472	c-----------------------------------------------------------------------
	473	c 3.4 Read Soil layers depths
	474	c-----------------------------------------------------------------------
	475
	476	ierr=NF_INQ_VARID(nid,"soildepth",nvarid)
	477	if (ierr.ne.NF_NOERR) then
	478	write(,)'lect_start_archive: Could not find <soildepth>'
	479	write(,)' => Assuming this is an archive in old format'
	480	olddepthdef=.true.
	481	depthinterpol=.true.
	482	! this is how soil depth was defined in ye old days
	483	do isoil=1,nsoilold
	484	mlayerold(isoil)=sqrt(887.75/3.14)((2.*(isoil-0.5))-1.)
	485	enddo
	486	else
	487	#ifdef NC_DOUBLE
	488	ierr = NF_GET_VAR_DOUBLE(nid,nvarid,mlayerold)
	489	#else
	490	ierr = NF_GET_VAR_REAL(nid,nvarid,mlayerold)
	491	#endif
	492	if (ierr .NE. NF_NOERR) then
	493	PRINT*, "lect_start_archive: Failed reading <soildepth>"
	494	CALL abort
	495	endif
	496
	497	endif !of if(ierr.ne.NF_NOERR)
	498
	499	! Read (or build) mlayer()
	500	if (depthinterpol) then
	501	! Build (default) new soil depths (mlayer(:) is in comsoil.h),
	502	! as in soil_settings.F
	503	write(,)' => Building default soil depths'
	504	do isoil=0,nsoilmx-1
	505	mlayer(isoil)=2.e-4(2.*(isoil-0.5))
	506	enddo
	507	write(,)' => mlayer: ',mlayer
	508	! Also build (default) new soil interlayer depth layer(:)
	509	do isoil=1,nsoilmx
	510	layer(isoil)=sqrt(mlayer(0)mlayer(1))
	511	& ((mlayer(1)/mlayer(0))**(isoil-1))
	512	enddo
	513	write(,)' => layer: ',layer
	514	else ! read mlayer() from file
	515	#ifdef NC_DOUBLE
	516	ierr = NF_GET_VAR_DOUBLE(nid,nvarid,mlayer)
	517	#else
	518	ierr = NF_GET_VAR_REAL(nid,nvarid,mlayer)
	519	#endif
	520	if (ierr .NE. NF_NOERR) then
	521	PRINT*, "lect_start_archive: Failed reading <soildepth>"
	522	CALL abort
	523	endif
[311]	524	! Also build (default) soil interlayer depth layer(:)
	525	do isoil=1,nsoilmx
	526	layer(isoil)=sqrt(mlayer(0)mlayer(1))
	527	& ((mlayer(1)/mlayer(0))**(isoil-1))
	528	enddo
[38]	529	endif ! of if (depthinterpol)
	530
	531	c-----------------------------------------------------------------------
	532	c 3.5 Read Soil thermal inertia
	533	c-----------------------------------------------------------------------
	534
	535	ierr=NF_INQ_VARID(nid,"inertiedat",nvarid)
	536	if (ierr.ne.NF_NOERR) then
	537	write(,)'lect_start_archive: Could not find <inertiedat>'
	538	write(,)' => Assuming this is an archive in old format'
	539	therminertia_3D=.false.
	540	write(,)' => Thermal inertia will be read from reference file'
	541	volcapa=1.e6
	542	write(,)' and soil volumetric heat capacity is set to ',
	543	& volcapa
	544	else
	545	#ifdef NC_DOUBLE
	546	ierr = NF_GET_VAR_DOUBLE(nid,nvarid,inertiedatold)
	547	#else
	548	ierr = NF_GET_VAR_REAL(nid,nvarid,inertiedatold)
	549	#endif
	550	if (ierr .NE. NF_NOERR) then
	551	PRINT*, "lect_start_archive: Failed reading <inertiedat>"
	552	CALL abort
	553	endif
	554	endif
	555
	556	c-----------------------------------------------------------------------
	557	c 3.6 Lecture geopotentiel au sol
	558	c-----------------------------------------------------------------------
	559	c
	560	ierr = NF_INQ_VARID (nid, "phisinit", nvarid)
	561	IF (ierr .NE. NF_NOERR) THEN
[146]	562	PRINT*, "lect_start_archive: <phisinit> is missing"
[38]	563	CALL abort
	564	ENDIF
	565	#ifdef NC_DOUBLE
	566	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, phisold)
	567	#else
	568	ierr = NF_GET_VAR_REAL(nid, nvarid, phisold)
	569	#endif
	570	IF (ierr .NE. NF_NOERR) THEN
[146]	571	PRINT*, "lect_start_archive: Failed loading <phisinit>"
[38]	572	CALL abort
	573	ENDIF
	574
	575	C-----------------------------------------------------------------------
	576	c lecture de "ptotalold" et "co2icetotalold"
	577	c-----------------------------------------------------------------------
	578	ptotalold = tab_cntrl(tab0+49)
	579	co2icetotalold = tab_cntrl(tab0+50)
	580
	581	c-----------------------------------------------------------------------
	582	c 4. Lecture du temps et choix
	583	c-----------------------------------------------------------------------
	584
	585	c lecture du temps
	586	c
	587	ierr = NF_INQ_DIMID (nid, "Time", nvarid)
	588	IF (ierr .NE. NF_NOERR) THEN
	589	ierr = NF_INQ_DIMID (nid, "temps", nvarid)
	590	IF (ierr .NE. NF_NOERR) THEN
[146]	591	PRINT*, "lect_start_archive: <Time> is missing"
[38]	592	CALL abort
	593	endif
	594	ENDIF
	595
	596	ierr = NF_INQ_VARID (nid, "Time", nvarid)
	597	IF (ierr .NE. NF_NOERR) THEN
	598	ierr = NF_INQ_VARID (nid, "temps", nvarid)
	599	endif
	600	#ifdef NC_DOUBLE
	601	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, timelist)
	602	#else
	603	ierr = NF_GET_VAR_REAL(nid, nvarid, timelist)
	604	#endif
	605	IF (ierr .NE. NF_NOERR) THEN
[146]	606	PRINT*, "lect_start_archive: Failed loading <Time>"
[38]	607	CALL abort
	608	ENDIF
	609	c
	610	write(,)
	611	write(,)
[146]	612	write(,) 'Dates of the stored initial states:'
	613	write(,) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
[38]	614	pi=2.*ASIN(1.)
	615	do i=1,timelen
	616	c call solarlong(timelist(i),sollong(i))
	617	c sollong(i) = sollong(i)*180./pi
[999]	618	c write(,) 'initial state at martian day: ',int(timelist(i))
	619	write(,) 'initial state at martian day: ',timelist(i)
[38]	620	c write(*,6) nint(timelist(i)),nint(mod(timelist(i),669)),
	621	c . sollong(i)
	622	end do
	623
	624	6 FORMAT(i7,i7,f9.3)
	625
	626	write(,)
[146]	627	write(,) 'Choose the martian day to use'
[38]	628	123 read(,,iostat=ierr) date
	629	if(ierr.ne.0) goto 123
	630	memo = 0
	631	do i=1,timelen
[999]	632	c if (date.eq.int(timelist(i))) then
	633	if (abs(date-timelist(i)).lt.0.01) then
[38]	634	memo = i
	635	endif
	636	end do
	637
	638	if (memo.eq.0) then
	639	write(,)
	640	write(,)
[146]	641	write(,) 'Wrong value for day number !!'
[38]	642	write(,)
[146]	643	write(,) 'Dates of the stored initial states:'
	644	write(,) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
[38]	645	do i=1,timelen
[999]	646	write(,) 'initial state at martian day: ',timelist(i)
[38]	647	c write(*,6) nint(timelist(i)),nint(mod(timelist(i),669))
	648	end do
	649	goto 123
	650	endif
[999]	651
[38]	652	!-----------------------------------------------------------------------
	653	! 5. Read (time-dependent) data from datafile
	654	!-----------------------------------------------------------------------
	655
	656
	657	c-----------------------------------------------------------------------
	658	c 5.1 Lecture des champs 2D (co2ice, emis,ps,tsurf,Tg[10], q2surf)
	659	c-----------------------------------------------------------------------
	660
	661	start=(/1,1,memo,0/)
	662	count=(/imold+1,jmold+1,1,0/)
	663
	664	ierr = NF_INQ_VARID (nid, "co2ice", nvarid)
	665	IF (ierr .NE. NF_NOERR) THEN
[146]	666	PRINT*, "lect_start_archive: <co2ice> is missing"
[38]	667	CALL abort
	668	ENDIF
	669	#ifdef NC_DOUBLE
	670	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,co2iceold)
	671	#else
	672	ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,co2iceold)
	673	#endif
	674	IF (ierr .NE. NF_NOERR) THEN
[146]	675	PRINT*, "lect_start_archive: Failed loading <co2ice>"
[38]	676	PRINT*, NF_STRERROR(ierr)
	677	CALL abort
	678	ENDIF
	679	c
	680	ierr = NF_INQ_VARID (nid, "emis", nvarid)
	681	IF (ierr .NE. NF_NOERR) THEN
[146]	682	PRINT*, "lect_start_archive: <emis> is missing"
[38]	683	CALL abort
	684	ENDIF
	685	#ifdef NC_DOUBLE
	686	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,emisold)
	687	#else
	688	ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,emisold)
	689	#endif
	690	IF (ierr .NE. NF_NOERR) THEN
[146]	691	PRINT*, "lect_start_archive: Failed loading <emis>"
[38]	692	CALL abort
	693	ENDIF
	694	c
	695	ierr = NF_INQ_VARID (nid, "ps", nvarid)
	696	IF (ierr .NE. NF_NOERR) THEN
[146]	697	PRINT*, "lect_start_archive: <ps> is missing"
[38]	698	CALL abort
	699	ENDIF
	700	#ifdef NC_DOUBLE
	701	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,psold)
	702	#else
	703	ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,psold)
	704	#endif
	705	IF (ierr .NE. NF_NOERR) THEN
[146]	706	PRINT*, "lect_start_archive: Failed loading <ps>"
[38]	707	CALL abort
	708	ENDIF
	709	c
	710	ierr = NF_INQ_VARID (nid, "tsurf", nvarid)
	711	IF (ierr .NE. NF_NOERR) THEN
[146]	712	PRINT*, "lect_start_archive: <tsurf> is missing"
[38]	713	CALL abort
	714	ENDIF
	715	#ifdef NC_DOUBLE
	716	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,tsurfold)
	717	#else
	718	ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,tsurfold)
	719	#endif
	720	IF (ierr .NE. NF_NOERR) THEN
[146]	721	PRINT*, "lect_start_archive: Failed loading <tsurf>"
[38]	722	CALL abort
	723	ENDIF
	724	c
	725	ierr = NF_INQ_VARID (nid, "q2surf", nvarid)
	726	IF (ierr .NE. NF_NOERR) THEN
[146]	727	PRINT*, "lect_start_archive: <q2surf> is missing"
[38]	728	CALL abort
	729	ENDIF
	730	#ifdef NC_DOUBLE
	731	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,q2old)
	732	#else
	733	ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,q2old)
	734	#endif
	735	IF (ierr .NE. NF_NOERR) THEN
[146]	736	PRINT*, "lect_start_archive: Failed loading <q2surf>"
[38]	737	CALL abort
	738	ENDIF
	739	c
	740	write(,)"lect_start_archive: rlonuold:"
	741	& ,rlonuold," rlatvold:",rlatvold
	742	write(,)
	743	c
	744
	745	c tracers: the 2 last ones are kept the 2 last one.
	746	c the others keep their rank. ! No longer true.
	747	c -------------------------------------------
	748	! Surface tracers:
[1036]	749	qsurfold(1:imold+1,1:jmold+1,1:nqtot)=0
[38]	750
[1036]	751	DO iq=1,nqtot
[38]	752	IF (oldtracernames) THEN
	753	txt=" "
	754	write(txt,'(a5,i2.2)')'qsurf',iq
	755	ELSE
	756	txt=trim(tnom(iq))//"_surf"
	757	if (txt.eq."h2o_vap") then
	758	! There is no surface tracer for h2o_vap;
	759	! "h2o_ice" should be loaded instead
	760	txt="h2o_ice_surf"
	761	write(,) 'lect_start_archive: loading surface tracer',
	762	& ' h2o_ice instead of h2o_vap'
	763	endif
	764	ENDIF ! of IF (oldtracernames)
	765	write(,) "lect_start_archive: loading tracer ",trim(txt)
	766	ierr = NF_INQ_VARID (nid,txt,nvarid)
	767	IF (ierr .NE. NF_NOERR) THEN
	768	PRINT*, "lect_start_archive: ",
	769	& " Tracer <",trim(txt),"> not found"
[146]	770	print*, "which (constant) value should it be initialized to?"
	771	read(,) tmpval
	772	qsurfold(1:imold+1,1:jmold+1,iq)=tmpval
[563]	773	ELSE ! tracer exists in file, load it
[38]	774	#ifdef NC_DOUBLE
[563]	775	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,
[38]	776	& qsurfold(1,1,iq))
	777	#else
[563]	778	ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,
[38]	779	& qsurfold(1,1,iq))
	780	#endif
[563]	781	IF (ierr .NE. NF_NOERR) THEN
	782	PRINT*, "lect_start_archive: ",
[38]	783	& " Failed loading <",trim(txt),">"
[563]	784	stop
	785	ENDIF
[38]	786	ENDIF
	787
[1036]	788	ENDDO ! of DO iq=1,nqtot
[38]	789
	790	!-----------------------------------------------------------------------
	791	! 5.2 Read 3D subterranean fields
	792	!-----------------------------------------------------------------------
	793
	794	start=(/1,1,1,memo/)
	795	count=(/imold+1,jmold+1,nsoilold,1/)
	796	!
	797	! Read soil temperatures
	798	!
	799	if (olddepthdef) then ! tsoil stored using the 'old format'
	800	start=(/1,1,memo,0/)
	801	count=(/imold+1,jmold+1,1,0/) ! because the "Tg" are 2D datasets
	802	do isoil=1,nsoilold
	803	! write(,)'isoil',isoil
	804	write(str2,'(i2.2)') isoil
	805	c
	806	ierr = NF_INQ_VARID (nid, "Tg"//str2, nvarid)
	807	IF (ierr .NE. NF_NOERR) THEN
	808	PRINT*, "lect_start_archive: ",
	809	& "Field <","Tg"//str2,"> not found"
	810	CALL abort
	811	ENDIF
	812	#ifdef NC_DOUBLE
	813	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,
	814	& tsoilold(1,1,isoil))
	815	#else
	816	ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,
	817	& tsoilold(1,1,isoil))
	818	#endif
	819	IF (ierr .NE. NF_NOERR) THEN
	820	PRINT*, "lect_start_archive: ",
	821	& "Failed reading <","Tg"//str2,">"
	822	CALL abort
	823	ENDIF
	824	c
	825	enddo ! of do isoil=1,nsoilold
	826
	827	! reset 'start' and 'count' to "3D" behaviour
	828	start=(/1,1,1,memo/)
	829	count=(/imold+1,jmold+1,nsoilold,1/)
	830
	831	else
	832	write(,) "lect_start_archive: loading tsoil "
	833	ierr=NF_INQ_VARID(nid,"tsoil",nvarid)
	834	if (ierr.ne.NF_NOERR) then
	835	write(,)"lect_start_archive: Cannot find <tsoil>"
	836	call abort
	837	else
	838	#ifdef NC_DOUBLE
	839	ierr=NF_GET_VARA_DOUBLE(nid,nvarid,start,count,tsoilold)
	840	#else
	841	ierr=NF_GET_VARA_REAL(nid,nvarid,start,count,tsoilold)
	842	#endif
	843	endif ! of if (ierr.ne.NF_NOERR)
	844
	845	endif ! of if (olddepthdef)
	846
	847	!
	848	! Read soil thermal inertias
	849	!
	850	! if (.not.olddepthdef) then ! no thermal inertia data in "old" archives
	851	! ierr=NF_INQ_VARID(nid,"inertiedat",nvarid)
	852	! if (ierr.ne.NF_NOERR) then
	853	! write(,)"lect_start_archive: Cannot find <inertiedat>"
	854	! call abort
	855	! else
	856	!#ifdef NC_DOUBLE
	857	! ierr=NF_GET_VARA_DOUBLE(nid,nvarid,start,count,inertiedatold)
	858	!#else
	859	! ierr=NF_GET_VARA_REAL(nid,nvarid,start,count,inertiedatold)
	860	!#endif
	861	! endif ! of if (ierr.ne.NF_NOERR)
	862	! endif
	863
	864	c-----------------------------------------------------------------------
	865	c 5.3 Lecture des champs 3D (t,u,v, q2atm,q)
	866	c-----------------------------------------------------------------------
	867
	868	start=(/1,1,1,memo/)
	869	count=(/imold+1,jmold+1,lmold,1/)
	870
	871	c
	872	ierr = NF_INQ_VARID (nid,"temp", nvarid)
	873	IF (ierr .NE. NF_NOERR) THEN
[146]	874	PRINT*, "lect_start_archive: <temp> is missing"
[38]	875	CALL abort
	876	ENDIF
	877	#ifdef NC_DOUBLE
	878	ierr = NF_GET_VARA_DOUBLE(nid, nvarid, start, count, told)
	879	#else
	880	ierr = NF_GET_VARA_REAL(nid, nvarid, start, count, told)
	881	#endif
	882	IF (ierr .NE. NF_NOERR) THEN
[146]	883	PRINT*, "lect_start_archive: Failed loading <temp>"
[38]	884	CALL abort
	885	ENDIF
	886	c
	887	ierr = NF_INQ_VARID (nid,"u", nvarid)
	888	IF (ierr .NE. NF_NOERR) THEN
[146]	889	PRINT*, "lect_start_archive: <u> is missing"
[38]	890	CALL abort
	891	ENDIF
	892	#ifdef NC_DOUBLE
	893	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,uold)
	894	#else
	895	ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,uold)
	896	#endif
	897	IF (ierr .NE. NF_NOERR) THEN
[146]	898	PRINT*, "lect_start_archive: Failed loading <u>"
[38]	899	CALL abort
	900	ENDIF
	901	c
	902	ierr = NF_INQ_VARID (nid,"v", nvarid)
	903	IF (ierr .NE. NF_NOERR) THEN
[146]	904	PRINT*, "lect_start_archive: <v> is missing"
[38]	905	CALL abort
	906	ENDIF
	907	#ifdef NC_DOUBLE
	908	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,vold)
	909	#else
	910	ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,vold)
	911	#endif
	912	IF (ierr .NE. NF_NOERR) THEN
[146]	913	PRINT*, "lect_start_archive: Failed loading <v>"
[38]	914	CALL abort
	915	ENDIF
	916	c
	917	ierr = NF_INQ_VARID (nid,"q2atm", nvarid)
	918	IF (ierr .NE. NF_NOERR) THEN
[146]	919	PRINT*, "lect_start_archive: <q2atm> is missing"
[38]	920	CALL abort
	921	ENDIF
	922	#ifdef NC_DOUBLE
	923	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,q2old(1,1,2))
	924	#else
	925	ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,q2old(1,1,2))
	926	#endif
	927	IF (ierr .NE. NF_NOERR) THEN
[146]	928	PRINT*, "lect_start_archive: Failed loading <q2atm>"
[38]	929	CALL abort
	930	ENDIF
	931	c
	932
	933	c tracers: the 2 last ones are kept the 2 last one.
	934	c the others keep their rank. ! No longer true.
	935	c -------------------------------------------
[146]	936	! Tracers:
[1036]	937	qold(1:imold+1,1:jmold+1,1:lmold,1:nqtot)=0
[38]	938
[1036]	939	DO iq=1,nqtot
[38]	940	IF (oldtracernames) THEN
	941	txt=" "
	942	write(txt,'(a1,i2.2)')'q',iq
	943	ELSE
	944	txt=tnom(iq)
	945	ENDIF
	946	write(,)"lect_start_archive: loading tracer ",trim(txt)
	947	ierr = NF_INQ_VARID (nid,txt,nvarid)
	948	IF (ierr .NE. NF_NOERR) THEN
	949	PRINT*, "lect_start_archive: ",
	950	& " Tracer <",trim(txt),"> not found"
[146]	951	print*, "which (constant) value should it be initialized to?"
	952	read(,) tmpval
	953	qold(1:imold+1,1:jmold+1,1:lmold,iq)=tmpval
[563]	954	ELSE ! tracer exists in file, load it
[38]	955	#ifdef NC_DOUBLE
[563]	956	ierr=NF_GET_VARA_DOUBLE(nid,nvarid,start,count,qold(1,1,1,iq))
[38]	957	#else
[563]	958	ierr=NF_GET_VARA_REAL(nid,nvarid,start,count,qold(1,1,1,iq))
[38]	959	#endif
[563]	960	IF (ierr .NE. NF_NOERR) THEN
	961	PRINT*, "lect_start_archive: ",
[38]	962	& " Failed loading <",trim(txt),">"
[563]	963	stop
	964	ENDIF
[38]	965	ENDIF
	966
[1036]	967	ENDDO ! of DO iq=1,nqtot
[38]	968
	969
	970	c Chemin pour trouver les donnees de surface (albedo, relief, th.inertia...)
	971	c -------------------------------------------------------------------------
	972
	973	! datapath = '/users/forget/gcm/data_mars_gcm'
	974
	975
	976	!=======================================================================
	977	! 6. Interpolation from old grid to new grid
	978	!=======================================================================
	979
	980	c=======================================================================
	981	c INTERPOLATION DANS LA NOUVELLE GRILLE et initialisation des variables
	982	c=======================================================================
	983	c Interpolation horizontale puis passage dans la grille physique pour
	984	c les variables physique
	985	c Interpolation verticale puis horizontale pour chaque variable 3D
	986	c=======================================================================
	987
	988	c-----------------------------------------------------------------------
	989	c 6.1 Variable 2d :
	990	c-----------------------------------------------------------------------
	991	c Relief
	992	call interp_horiz (phisold,phisold_newgrid,imold,jmold,iim,jjm,1,
	993	& rlonuold,rlatvold,rlonu,rlatv)
	994
	995	c Glace CO2
	996	call interp_horiz (co2iceold,co2ices,imold,jmold,iim,jjm,1,
	997	& rlonuold,rlatvold,rlonu,rlatv)
	998
	999	c Temperature de surface
	1000	call interp_horiz (tsurfold,tsurfs,imold,jmold,iim,jjm,1,
	1001	& rlonuold,rlatvold,rlonu,rlatv)
	1002	call gr_dyn_fi (1,iim+1,jjm+1,ngridmx,tsurfs,tsurf)
	1003	c write(44,*) 'tsurf', tsurf
	1004
	1005	c Temperature du sous-sol
	1006	! call interp_horiz(tsoilold,tsoils,
	1007	! & imold,jmold,iim,jjm,nsoilmx,
	1008	! & rlonuold,rlatvold,rlonu,rlatv)
	1009	! call gr_dyn_fi (nsoilmx,iim+1,jjm+1,ngridmx,tsoils,tsoil)
	1010	c write(45,*) 'tsoil',tsoil
	1011
	1012	c Emissivite de la surface
	1013	call interp_horiz (emisold,emiss,imold,jmold,iim,jjm,1,
	1014	& rlonuold,rlatvold,rlonu,rlatv)
	1015	call gr_dyn_fi (1,iim+1,jjm+1,ngridmx,emiss,emis)
	1016	c write(46,*) 'emis',emis
	1017	c-----------------------------------------------------------------------
	1018	c 6.1.2 Traitement special de la pression au sol :
	1019	c-----------------------------------------------------------------------
	1020
	1021	c Extrapolation la pression dans la nouvelle grille
	1022	call interp_horiz(psold,ps,imold,jmold,iim,jjm,1,
	1023	& rlonuold,rlatvold,rlonu,rlatv)
	1024
	1025	c-----------------------------------------------------------------------
	1026	c On assure la conservation de la masse de l'atmosphere + calottes
	1027	c-----------------------------------------------------------------------
	1028
	1029	ptotal = 0.
	1030	co2icetotal = 0.
	1031	DO j=1,jjp1
	1032	DO i=1,iim
	1033	ptotal=ptotal+ps(i,j)*aire(i,j)/g
	1034	co2icetotal = co2icetotal + co2iceS(i,j)*aire(i,j)
	1035	END DO
	1036	END DO
	1037
	1038	write(,)
[146]	1039	write(,)'Old grid: mass of the atmosphere :',ptotalold
	1040	write(,)'New grid: mass of the atmosphere :',ptotal
	1041	write (,) 'Ratio new atm / old atm =', ptotal/ptotalold
[38]	1042	write(,)
[146]	1043	write(,)'Old grid: mass of CO2 ice:',co2icetotalold
	1044	write(,)'New grid: mass of CO2 ice:',co2icetotal
[563]	1045	if (co2icetotalold.ne.0.) then
	1046	write(,)'Ratio new ice / old ice =',co2icetotal/co2icetotalold
	1047	endif
[38]	1048	write(,)
	1049
	1050
	1051	DO j=1,jjp1
	1052	DO i=1,iip1
	1053	ps(i,j)=ps(i,j) * ptotalold/ptotal
	1054	END DO
	1055	END DO
	1056
	1057	if ( co2icetotalold.gt.0.) then
	1058	DO j=1,jjp1
	1059	DO i=1,iip1
	1060	co2iceS(i,j)=co2iceS(i,j) * co2icetotalold/co2icetotal
	1061	END DO
	1062	END DO
	1063	end if
	1064
	1065	c-----------------------------------------------------------------------
	1066	c 6.2 Subterranean 3d variables:
	1067	c-----------------------------------------------------------------------
	1068
	1069	c-----------------------------------------------------------------------
	1070	c 6.2.1 Thermal Inertia
	1071	c Note: recall that inertiedat is a common in "comsoil.h"
	1072	c-----------------------------------------------------------------------
	1073
	1074	! depth-wise interpolation, if required
	1075	if (depthinterpol.and.(.not.olddepthdef)) then
	1076	allocate(oldval(nsoilold))
	1077	allocate(newval(nsoilmx))
	1078	write(,)'lect_start_archive: WARNING: vertical interpolation o
	1079	&f soil thermal inertia; might be wiser to reset it.'
	1080	write(,)
	1081
	1082	do i=1,imold+1
	1083	do j=1,jmold+1
	1084	!copy old values
	1085	oldval(1:nsoilold)=inertiedatold(i,j,1:nsoilold)
	1086	!interpolate
	1087	call interp_line(mlayerold,oldval,nsoilold,
	1088	& mlayer,newval,nsoilmx)
	1089	!copy interpolated values
	1090	inertiedatoldnew(i,j,1:nsoilmx)=newval(1:nsoilmx)
	1091	enddo
	1092	enddo
	1093	! cleanup
	1094	deallocate(oldval)
	1095	deallocate(newval)
	1096	endif !of if (depthinterpol)
	1097
	1098	if (therminertia_3D) then
	1099	! We have inertiedatold
	1100	if((imold.ne.iim).or.(jmold.ne.jjm)) then
	1101	write(,)'lect_start_archive: WARNING: horizontal interpolation
	1102	&of thermal inertia; might be better to reset it.'
	1103	write(,)
	1104	endif
	1105
	1106	! Do horizontal interpolation
	1107	if (depthinterpol) then
	1108	call interp_horiz(inertiedatoldnew,inertiedatS,
	1109	& imold,jmold,iim,jjm,nsoilmx,
	1110	& rlonuold,rlatvold,rlonu,rlatv)
	1111	else
	1112	call interp_horiz(inertiedatold,inertiedatS,
	1113	& imold,jmold,iim,jjm,nsoilold,
	1114	& rlonuold,rlatvold,rlonu,rlatv)
	1115	endif ! of if (depthinterpol)
	1116
	1117	else ! no 3D thermal inertia data
	1118	write(,)'lect_start_archive: using reference surface inertia'
	1119	! Use surface inertia (and extend it to all depths)
	1120	do i=1,nsoilmx
	1121	inertiedatS(1:iip1,1:jjp1,i)=surfith(1:iip1,1:jjp1)
	1122	enddo
	1123	! Build an old resolution surface thermal inertia
	1124	! (will be needed for tsoil interpolation)
	1125	call interp_horiz(surfith,surfithold,
	1126	& iim,jjm,imold,jmold,1,
	1127	& rlonu,rlatv,rlonuold,rlatvold)
	1128	endif
	1129
	1130
	1131	! Reshape inertiedatS to scalar grid as inertiedat
	1132	call gr_dyn_fi (nsoilmx,iim+1,jjm+1,ngridmx,
	1133	& inertiedatS,inertiedat)
	1134
	1135	c-----------------------------------------------------------------------
	1136	c 6.2.2 Soil temperature
	1137	c-----------------------------------------------------------------------
	1138	! write(,) 'Soil'
	1139	! Recast temperatures along soil depth, if necessary
	1140	if (olddepthdef) then
	1141	allocate(oldgrid(nsoilold+1))
	1142	allocate(oldval(nsoilold+1))
	1143	allocate(newval(nsoilmx))
[174]	1144	do i=1,imold+1
	1145	do j=1,jmold+1
[38]	1146	! copy values
[146]	1147	oldval(1)=tsurfold(i,j)
[38]	1148	oldval(2:nsoilold+1)=tsoilold(i,j,1:nsoilold)
	1149	! build vertical coordinate
	1150	oldgrid(1)=0. ! ground
	1151	oldgrid(2:nsoilold+1)=mlayerold(1:nsoilold)*
	1152	& (surfithold(i,j)/1.e6)
	1153	! Note; at this stage, we impose volcapa=1.e6 above
	1154	! since volcapa isn't set in old soil definitions
	1155
	1156	! interpolate
	1157	call interp_line(oldgrid,oldval,nsoilold+1,
	1158	& mlayer,newval,nsoilmx)
	1159	! copy result in tsoilold
	1160	tsoiloldnew(i,j,1:nsoilmx)=newval(1:nsoilmx)
	1161	enddo
	1162	enddo
	1163	! cleanup
	1164	deallocate(oldgrid)
	1165	deallocate(oldval)
	1166	deallocate(newval)
	1167
	1168	else
	1169	if (depthinterpol) then ! if vertical interpolation is required
	1170	allocate(oldgrid(nsoilold+1))
	1171	allocate(oldval(nsoilold+1))
	1172	allocate(newval(nsoilmx))
	1173	! build vertical coordinate
	1174	oldgrid(1)=0. ! ground
	1175	oldgrid(2:nsoilold+1)=mlayerold(1:nsoilold)
[174]	1176	do i=1,imold+1
	1177	do j=1,jmold+1
[38]	1178	! copy values
[146]	1179	oldval(1)=tsurfold(i,j)
[38]	1180	oldval(2:nsoilold+1)=tsoilold(i,j,1:nsoilold)
	1181	! interpolate
	1182	call interp_line(oldgrid,oldval,nsoilold+1,
	1183	& mlayer,newval,nsoilmx)
	1184	! copy result in tsoilold
	1185	tsoiloldnew(i,j,1:nsoilmx)=newval(1:nsoilmx)
	1186	enddo
	1187	enddo
	1188	! write(,)'tsoiloldnew(1,1,1):',tsoiloldnew(1,1,1)
	1189	! cleanup
	1190	deallocate(oldgrid)
	1191	deallocate(oldval)
	1192	deallocate(newval)
	1193
	1194	else
	1195	tsoiloldnew(:,:,:)=tsoilold(:,:,:)
	1196	endif ! of if (depthinterpol)
	1197	endif ! of if (olddepthdef)
	1198
	1199	! Do the horizontal interpolation
	1200	call interp_horiz(tsoiloldnew,tsoilS,
	1201	& imold,jmold,iim,jjm,nsoilmx,
	1202	& rlonuold,rlatvold,rlonu,rlatv)
	1203
	1204	! Reshape tsoilS to scalar grid as tsoil
	1205	call gr_dyn_fi (nsoilmx,iim+1,jjm+1,ngridmx,tsoilS,tsoil)
	1206
	1207
	1208
	1209	c-----------------------------------------------------------------------
	1210	c 6.3 Variable 3d :
	1211	c-----------------------------------------------------------------------
	1212
	1213	c temperatures atmospheriques
	1214	write (,) 'lect_start_archive: told ', told (1,jmold+1,1) ! INFO
	1215	call interp_vert
	1216	& (told,var,lmold,llm,apsold,bpsold,aps,bps,
	1217	& psold,(imold+1)*(jmold+1))
	1218	write (,) 'lect_start_archive: var ', var (1,jmold+1,1) ! INFO
	1219	call interp_horiz(var,t,imold,jmold,iim,jjm,llm,
	1220	& rlonuold,rlatvold,rlonu,rlatv)
	1221	write (,) 'lect_start_archive: t ', t(1,jjp1,1) ! INFO
	1222
	1223	c q2 : pbl wind variance
	1224	write (,) 'lect_start_archive: q2old ', q2old (1,2,1) ! INFO
	1225	call interp_vert (q2old,varp1,lmold+1,llm+1,
	1226	& apsold,bpsold,aps,bps,psold,(imold+1)*(jmold+1))
	1227	write (,) 'lect_start_archive: varp1 ', varp1 (1,2,1) ! INFO
	1228	call interp_horiz(varp1,q2s,imold,jmold,iim,jjm,llm+1,
	1229	& rlonuold,rlatvold,rlonu,rlatv)
	1230	write (,) 'lect_start_archive: q2s ', q2s (1,2,1) ! INFO
	1231	call gr_dyn_fi (llm+1,iim+1,jjm+1,ngridmx,q2s,q2)
	1232	write (,) 'lect_start_archive: q2 ', q2 (1,2) ! INFO
	1233	c write(47,*) 'q2',q2
	1234
	1235	c calcul des champ de vent; passage en vent covariant
	1236	write (,) 'lect_start_archive: uold ', uold (1,2,1) ! INFO
	1237	call interp_vert
	1238	& (uold,var,lmold,llm,apsold,bpsold,aps,bps,
	1239	& psold,(imold+1)*(jmold+1))
	1240	write (,) 'lect_start_archive: var ', var (1,2,1) ! INFO
	1241	call interp_horiz(var,us,imold,jmold,iim,jjm,llm,
	1242	& rlonuold,rlatvold,rlonu,rlatv)
	1243	write (,) 'lect_start_archive: us ', us (1,2,1) ! INFO
	1244
	1245	call interp_vert
	1246	& (vold,var,lmold,llm,
	1247	& apsold,bpsold,aps,bps,psold,(imold+1)*(jmold+1))
	1248	call interp_horiz(var,vs,imold,jmold,iim,jjm,llm,
	1249	& rlonuold,rlatvold,rlonu,rlatv)
	1250	call scal_wind(us,vs,unat,vnat)
	1251	write (,) 'lect_start_archive: unat ', unat (1,2,1) ! INFO
	1252	do l=1,llm
	1253	do j = 1, jjp1
	1254	do i=1,iip1
	1255	ucov( i,j,l ) = unat( i,j,l ) * cu(i,j)
	1256	c ucov( i,j,l ) = 0
	1257	end do
	1258	end do
	1259	end do
	1260	write (,) 'lect_start_archive: ucov ', ucov (1,2,1) ! INFO
	1261	c write(48,*) 'ucov',ucov
	1262	do l=1,llm
	1263	do j = 1, jjm
	1264	do i=1,iim
	1265	vcov( i,j,l ) = vnat( i,j,l ) * cv(i,j)
	1266	c vcov( i,j,l ) = 0
	1267	end do
	1268	vcov( iip1,j,l ) = vcov( 1,j,l )
	1269	end do
	1270	end do
	1271	c write(49,*) 'ucov',vcov
	1272
	1273	c traceurs surface
[1036]	1274	do iq = 1, nqtot
[38]	1275	call interp_horiz(qsurfold(1,1,iq) ,qsurfs(1,1,iq),
	1276	& imold,jmold,iim,jjm,1,
	1277	& rlonuold,rlatvold,rlonu,rlatv)
	1278	enddo
	1279
[1036]	1280	call gr_dyn_fi (nqtot,iim+1,jjm+1,ngridmx,qsurfs,qsurf)
[38]	1281
	1282	c traceurs 3D
[1036]	1283	do iq = 1, nqtot
[38]	1284	call interp_vert(qold(1,1,1,iq),var,lmold,llm,
	1285	& apsold,bpsold,aps,bps,psold,(imold+1)*(jmold+1))
	1286	call interp_horiz(var,q(1,1,1,iq),imold,jmold,iim,jjm,llm,
	1287	& rlonuold,rlatvold,rlonu,rlatv)
	1288	enddo
	1289	cccccccccccccccccccccccccccccc
	1290	c make sure that sum of q = 1
	1291	c dominent species is = 1 - sum(all other species)
	1292	cccccccccccccccccccccccccccccc
	1293	c iqmax=1
	1294	c
	1295	c if (nqold.gt.10) then
	1296	c do l=1,llm
	1297	c do j=1,jjp1
	1298	c do i=1,iip1
	1299	c do iq=1,nqold
	1300	c if (q(i,j,l,iq).gt.q(i,j,l,iqmax)) then
	1301	c iqmax=iq
	1302	c endif
	1303	c enddo
	1304	c q(i,j,l,iqmax)=1.
	1305	c qtot(i,j,l)=0
	1306	c do iq=1,nqold
	1307	c if (iq.ne.iqmax) then
	1308	c q(i,j,l,iqmax)=q(i,j,l,iqmax)-q(i,j,l,iq)
	1309	c endif
	1310	c enddo !iq
	1311	c do iq=1,nqold
	1312	c qtot(i,j,l)=qtot(i,j,l)+q(i,j,l,iq)
	1313	c if (i.eq.1.and.j.eq.1.and.l.Eq.1) write(,)' qtot(i,j,l)',
	1314	c $ qtot(i,j,l)
	1315	c enddo !iq
	1316	c enddo !i
	1317	c enddo !j
	1318	c enddo !l
	1319	c endif
	1320	ccccccccccccccccccccccccccccccc
	1321
	1322	c Periodicite :
[1036]	1323	do iq = 1, nqtot
[38]	1324	do l=1, llm
	1325	do j = 1, jjp1
	1326	q(iip1,j,l,iq) = q(1,j,l,iq)
	1327	end do
	1328	end do
	1329	enddo
	1330
	1331	call gr_dyn_fi (1,iim+1,jjm+1,ngridmx,co2ices,co2ice)
	1332
	1333	c-----------------------------------------------------------------------
	1334	c Initialisation h: (passage de t -> h)
	1335	c-----------------------------------------------------------------------
	1336
	1337	DO l=1,llm
	1338	DO j=1,jjp1
	1339	DO i=1,iim
	1340	h(i,j,l) = t(i,j,l)((ps(i,j)/preff)*kappa)
	1341	END DO
	1342	h(iip1,j,l) = h(1,j,l)
	1343	END DO
	1344	END DO
	1345
	1346
	1347	c***********************************************************************
	1348	c***********************************************************************
	1349	c Fin subroutine lecture ini
	1350	c***********************************************************************
	1351	c***********************************************************************
	1352
	1353	deallocate(timelist)
	1354	deallocate(rlonuold, rlatvold)
	1355	deallocate(rlonvold, rlatuold)
	1356	deallocate(apsold,bpsold)
	1357	deallocate(uold)
	1358	deallocate(vold)
	1359	deallocate(told)
	1360	deallocate(psold)
	1361	deallocate(phisold)
	1362	deallocate(qold)
	1363	deallocate(co2iceold)
	1364	deallocate(tsurfold)
	1365	deallocate(emisold)
	1366	deallocate(q2old)
	1367	deallocate(tsoilold)
	1368	deallocate(tsoiloldnew)
	1369	deallocate(inertiedatold)
	1370	deallocate(inertiedatoldnew)
	1371	deallocate(surfithold)
	1372	deallocate(mlayerold)
	1373	deallocate(qsurfold)
	1374	deallocate(var,varp1)
	1375
	1376	! write(,)'lect_start_archive: END'
	1377	return
	1378	end

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

Download in other formats:

Original Format