Context Navigation

lect_start_archive.F @ 1202

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

Generic GCM:

Some more cleanup in dynamics:
- Moved "start2archive" (and auxilliary routines) to phystd
- removed unused (obsolete) testharm.F , para_netcdf.h , readhead_NC.F , angtot.h from dyn3d
- removed obsolete addit.F (and change corresponding lines in gcm)
- remove unused "description.h" (and many places where it was "included")

File size: 44.4 KB

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

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

Download in other formats:

Original Format