Context Navigation

lect_start_archive.F @ 1390

Last change on this file since 1390 was 1266, checked in by aslmd, 11 years ago

LMDZ.MARS
IMPORTANT CHANGE

Remove all reference/use of nlayermx and dimphys.h
Made use of automatic arrays whenever arrays are needed with dimension nlayer
Remove lots of obsolete reference to dimensions.h
Converted iono.h and param_v4.h into corresponding modules

(with embedded subroutine to allocate arrays)
(no arrays allocated if thermosphere not used)

Deleted param.h and put contents into module param_v4_h
Adapted testphys1d, newstart, etc...
Made DATA arrays in param_read to be initialized by subroutine

fill_data_thermos in module param_v4_h

Optimized computations in paramfoto_compact (twice less dlog10 calculations)
Checked consistency before/after modification in debug mode
Checked performance is not impacted (same as before)

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

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

Download in other formats:

Original Format