Context Navigation

lect_start_archive.F @ 621

Last change on this file since 621 was 563, checked in by emillour, 13 years ago

Mars GCM:

Added option "q=profile" in newstart to initialize a given tracer with

a profile specified in an ASCII file.

Minor fix on lect_start_archive : when tracers were not found, a default

value was asked to user, but twice! Once is enough.

File size: 44.8 KB

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

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

Download in other formats:

Original Format