Context Navigation

lect_start_archive.F @ 1212

Last change on this file since 1212 was 1208, checked in by tnavarro, 11 years ago
added tauscaling in startfi + moved start_archive routines from dyn3d to phymars
File size: 46.2 KB

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

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

Download in other formats:

Original Format