Context Navigation

newstart.F @ 3094

Last change on this file since 3094 was 3025, checked in by emillour, 17 months ago
Mars PCM: Fixed a bug in newstart call to lect_start_archive; missing perenial_co2ice argument. Turned lect_start_archive into a module so this cant't happen again. EM
File size: 76.1 KB

Rev	Line
[38]	1	C======================================================================
	2	PROGRAM newstart
	3	c=======================================================================
	4	c
	5	c
	6	c Auteur: Christophe Hourdin/Francois Forget/Yann Wanherdrick
	7	c ------
	8	c Derniere modif : 12/03
	9	c
	10	c
	11	c Objet: Create or modify the initial state for the LMD Mars GCM
	12	c ----- (fichiers NetCDF start et startfi)
	13	c
	14	c
	15	c=======================================================================
	16
[1543]	17	use ioipsl_getincom, only: getin
	18	use mod_phys_lmdz_para, only: is_parallel, is_sequential,
	19	& is_mpi_root, is_omp_root,
	20	& is_master
[1415]	21	use infotrac, only: infotrac_init, nqtot, tname
[1232]	22	use tracer_mod, only: noms, mmol,
	23	& igcm_dust_number, igcm_dust_mass,
[1130]	24	& igcm_ccn_number, igcm_ccn_mass,
[1232]	25	& igcm_h2o_vap, igcm_h2o_ice, igcm_co2,
[2312]	26	& igcm_hdo_vap, igcm_hdo_ice,
[2258]	27	& igcm_n2, igcm_ar, igcm_o2, igcm_co,
	28	& igcm_o, igcm_h2
[1047]	29	use surfdat_h, only: phisfi, z0, zmea, zstd, zsig, zgam, zthe,
[1232]	30	& albedodat, z0_default, qsurf, tsurf,
[2909]	31	& emis, hmons, summit, base, watercap,
[2999]	32	& ini_surfdat_h_slope_var,end_surfdat_h_slope_var,
	33	& perenial_co2ice
[2942]	34	use comsoil_h, only: inertiedat, inertiesoil,layer, mlayer,
[2952]	35	& nsoilmx,tsoil,ini_comsoil_h_slope_var, end_comsoil_h_slope_var,
	36	& flux_geo
[1415]	37	use control_mod, only: day_step, iphysiq, anneeref, planet_type
[2001]	38	use geometry_mod, only: longitude,latitude,cell_area
[3025]	39	use lect_start_archive_mod, only: lect_start_archive
[1944]	40	use phyetat0_mod, only: phyetat0
[1130]	41	use phyredem, only: physdem0, physdem1
	42	use iostart, only: open_startphy
[2909]	43	use dimradmars_mod, only: albedo,
	44	& ini_dimradmars_mod_slope_var,end_dimradmars_mod_slope_var
[2410]	45	use dust_param_mod, only: tauscaling
[1944]	46	use turb_mod, only: q2, wstar
[1403]	47	use filtreg_mod, only: inifilr
[1543]	48	USE mod_const_mpi, ONLY: COMM_LMDZ
[1422]	49	USE comvert_mod, ONLY: ap,bp,pa,preff
	50	USE comconst_mod, ONLY: lllm,daysec,dtphys,dtvr,
	51	. cpp,kappa,rad,omeg,g,r,pi
	52	USE serre_mod, ONLY: alphax
	53	USE temps_mod, ONLY: day_ini,hour_ini
	54	USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
[1543]	55	USE iniphysiq_mod, ONLY: iniphysiq
[2167]	56	USE exner_hyb_m, ONLY: exner_hyb
[2673]	57	USE inichim_newstart_mod, ONLY: inichim_newstart
[2909]	58	use comslope_mod, ONLY: nslope,def_slope,def_slope_mean,
	59	& subslope_dist,end_comslope_h,ini_comslope_h
[2910]	60	USE subslope_mola_mod, ONLY: subslope_mola
	61
[38]	62	implicit none
	63
[1543]	64	include "dimensions.h"
[1047]	65	integer, parameter :: ngridmx = (2+(jjm-1)*iim - 1/jjm)
[1543]	66	include "paramet.h"
	67	include "comgeom2.h"
	68	include "comdissnew.h"
	69	include "clesph0.h"
	70	include "netcdf.inc"
[38]	71	c=======================================================================
	72	c Declarations
	73	c=======================================================================
	74
	75	c Variables dimension du fichier "start_archive"
	76	c------------------------------------
	77	CHARACTER relief*3
	78
	79	c et autres:
	80	c----------
	81
	82	c Variables pour les lectures NetCDF des fichiers "start_archive"
	83	c--------------------------------------------------
	84	INTEGER nid_dyn, nid_fi,nid,nvarid
	85	INTEGER tab0
	86
	87	REAL date
	88	REAL p_rad,p_omeg,p_g,p_mugaz,p_daysec
	89
	90	c Variable histoire
	91	c------------------
	92	REAL vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm) ! vents covariants
	93	REAL phis(iip1,jjp1)
[1036]	94	REAL,ALLOCATABLE :: q(:,:,:,:) ! champs advectes
[38]	95
	96	c autre variables dynamique nouvelle grille
	97	c------------------------------------------
	98	REAL pks(iip1,jjp1)
	99	REAL w(iip1,jjp1,llm+1)
	100	REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
	101	! REAL dv(ip1jm,llm),du(ip1jmp1,llm)
	102	! REAL dh(ip1jmp1,llm),dp(ip1jmp1)
	103	REAL phi(iip1,jjp1,llm)
	104
	105	integer klatdat,klongdat
	106	PARAMETER (klatdat=180,klongdat=360)
	107
	108	c Physique sur grille scalaire
	109	c----------------------------
	110	real zmeaS(iip1,jjp1),zstdS(iip1,jjp1)
	111	real zsigS(iip1,jjp1),zgamS(iip1,jjp1),ztheS(iip1,jjp1)
[1974]	112	real hmonsS(iip1,jjp1)
	113	real summitS(iip1,jjp1)
[2079]	114	real baseS(iip1,jjp1)
[1974]	115	real zavgS(iip1,jjp1)
[224]	116	real z0S(iip1,jjp1)
[38]	117
	118	c variable physique
	119	c------------------
[1208]	120	REAL tauscadyn(iip1,jjp1) ! dust conversion factor on the dynamics grid
[38]	121	real alb(iip1,jjp1),albfi(ngridmx) ! albedos
	122	real ith(iip1,jjp1,nsoilmx),ithfi(ngridmx,nsoilmx) ! thermal inertia (3D)
	123	real surfith(iip1,jjp1),surfithfi(ngridmx) ! surface thermal inertia (2D)
[2001]	124	! REAL latfi(ngridmx),lonfi(ngridmx),airefi(ngridmx)
[38]	125
	126	INTEGER i,j,l,isoil,ig,idum
	127	real mugaz ! molar mass of the atmosphere
	128
	129	integer ierr !, nbetat
	130
	131	c Variables on the new grid along scalar points
	132	c------------------------------------------------------
	133	! REAL p(iip1,jjp1)
	134	REAL t(iip1,jjp1,llm)
	135	real phisold_newgrid(iip1,jjp1)
	136	REAL :: teta(iip1, jjp1, llm)
	137	REAL :: pk(iip1,jjp1,llm)
	138	REAL :: pkf(iip1,jjp1,llm)
	139	REAL :: ps(iip1, jjp1)
	140	REAL :: masse(iip1,jjp1,llm)
	141	REAL :: xpn,xps,xppn(iim),xpps(iim)
	142	REAL :: p3d(iip1, jjp1, llm+1)
	143	! REAL dteta(ip1jmp1,llm)
	144
	145	c Variable de l'ancienne grille
	146	c------------------------------
	147	real time
	148	real tab_cntrl(100)
	149	real tab_cntrl_bis(100)
	150
	151	c variables diverses
	152	c-------------------
[1130]	153	real choix_1 ! ==0 : read start_archive file ; ==1: read start files
[38]	154	character*80 fichnom
[618]	155	integer Lmodif,iq
	156	integer flagthermo, flagh2o
[38]	157	character modif*20
	158	real tsud,albsud,alb_bb,ith_bb,Tiso
	159	real ptoto,pcap,patm,airetot,ptotn,patmn
	160	! real ssum
	161	character*1 yes
	162	logical :: flagiso=.false. , flagps0=.false.
	163	real val, val2, val3 ! to store temporary variables
	164	real :: iceith=2000 ! thermal inertia of subterranean ice
	165	real :: iceithN,iceithS ! values of thermal inertias in N & S hemispheres
	166	integer iref,jref
	167
	168	INTEGER :: itau
[2312]	169	real DoverH !D/H ratio
[38]	170
[1231]	171	INTEGER :: numvanle
[563]	172	character(len=50) :: txt ! to store some text
[38]	173	integer :: count
[563]	174	real :: profile(llm+1) ! to store an atmospheric profile + surface value
[38]	175
	176	! MONS data:
	177	real :: MONS_Hdn(iip1,jjp1) ! Hdn: %WEH=Mass fraction of H2O
	178	real :: MONS_d21(iip1,jjp1) ! ice table "depth" (in kg/m2)
	179	! coefficient to apply to convert d21 to 'true' depth (m)
	180	real :: MONS_coeff
	181	real :: MONS_coeffS ! coeff for southern hemisphere
	182	real :: MONS_coeffN ! coeff for northern hemisphere
	183	! real,parameter :: icedepthmin=1.e-3 ! Ice begins at most at that depth
[1232]	184	! Reference position for composition
	185	real :: latref,lonref,dlatmin,dlonmin
	186	! Variable used to change composition
[2175]	187	real :: Svmr,Smmr,Smmr_old,Smmr_new,n,Sn
[1232]	188	real :: Mair_old,Mair_new,vmr_old,vmr_new
	189	real,allocatable :: coefvmr(:) ! Correction coefficient when changing composition
[2258]	190	real :: maxq
[2909]	191	integer :: iloc(1), iqmax, islope
[1711]	192	! sub-grid cloud fraction
	193	real :: totcloudfrac(ngridmx)
[2909]	194	!Variable to change the number of subslope
	195	REAL, ALLOCATABLE :: default_def_slope(:)
	196	REAL,ALLOCATABLE :: tsurf_old_slope(:,:) ! Surface temperature (K)
	197	REAL,ALLOCATABLE :: emis_old_slope(:,:) ! Thermal IR surface emissivity
	198	REAL,ALLOCATABLE :: qsurf_old_slope(:,:,:) ! tracer on surface (e.g. kg.m-2)
	199	REAL,ALLOCATABLE :: watercap_old_slope(:,:) ! Surface water ice (kg.m-2)
[2999]	200	REAL,ALLOCATABLE :: perenial_co2_old_slope(:,:) ! Surface water ice (kg.m-2)
[2909]	201	REAL,ALLOCATABLE :: tsoil_old_slope(:,:,:)
[2952]	202	REAL,ALLOCATABLE :: inertiesoil_old_slope(:,:,:)
[2909]	203	REAL,ALLOCATABLE :: albedo_old_slope(:,:,:) ! Surface albedo in each solar band
[2952]	204	REAL,ALLOCATABLE :: flux_geo_old_slope(:,:)
[2909]	205	integer :: iflat
	206	integer :: nslope_old, nslope_new
[38]	207
[2909]	208
	209
[38]	210	c sortie visu pour les champs dynamiques
	211	c---------------------------------------
	212	! INTEGER :: visuid
	213	! real :: time_step,t_ops,t_wrt
	214	! CHARACTER*80 :: visu_file
	215
	216	cpp = 744.499 ! for Mars, instead of 1004.70885 (Earth)
	217	preff = 610. ! for Mars, instead of 101325. (Earth)
	218	pa= 20 ! for Mars, instead of 500 (Earth)
[1415]	219	planet_type="mars"
[38]	220
[1543]	221	! initialize "serial/parallel" related stuff:
	222	! (required because we call tabfi() below, before calling iniphysiq)
	223	is_sequential=.true.
	224	is_parallel=.false.
	225	is_mpi_root=.true.
	226	is_omp_root=.true.
	227	is_master=.true.
	228
[1036]	229	! Load tracer number and names:
[1415]	230	call infotrac_init
[1036]	231	! allocate arrays
	232	allocate(q(iip1,jjp1,llm,nqtot))
[1232]	233	allocate(coefvmr(nqtot))
	234
[38]	235	c=======================================================================
	236	c Choice of the start file(s) to use
	237	c=======================================================================
	238
	239	write(,) 'From which kind of files do you want to create new',
	240	. 'start and startfi files'
	241	write(,) ' 0 - from a file start_archive'
	242	write(,) ' 1 - from files start and startfi'
	243
	244	c-----------------------------------------------------------------------
	245	c Open file(s) to modify (start or start_archive)
	246	c-----------------------------------------------------------------------
	247
	248	DO
	249	read(,,iostat=ierr) choix_1
	250	if ((choix_1 /= 0).OR.(choix_1 /=1)) EXIT
	251	ENDDO
	252
	253	c Open start_archive
	254	c ~~~~~~~~~~~~~~~~~~~~~~~~~~
	255	if (choix_1.eq.0) then
	256
	257	write(,) 'Creating start files from:'
	258	write(,) './start_archive.nc'
	259	write(,)
	260	fichnom = 'start_archive.nc'
	261	ierr = NF_OPEN (fichnom, NF_NOWRITE,nid)
	262	IF (ierr.NE.NF_NOERR) THEN
	263	write(6,*)' Problem opening file:',fichnom
	264	write(6,*)' ierr = ', ierr
	265	CALL ABORT
	266	ENDIF
	267	tab0 = 50
	268	Lmodif = 1
	269
	270	c OR open start and startfi files
	271	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	272	else
	273	write(,) 'Creating start files from:'
	274	write(,) './start.nc and ./startfi.nc'
	275	write(,)
	276	fichnom = 'start.nc'
	277	ierr = NF_OPEN (fichnom, NF_NOWRITE,nid_dyn)
	278	IF (ierr.NE.NF_NOERR) THEN
	279	write(6,*)' Problem opening file:',fichnom
	280	write(6,*)' ierr = ', ierr
	281	CALL ABORT
	282	ENDIF
	283
	284	fichnom = 'startfi.nc'
	285	ierr = NF_OPEN (fichnom, NF_NOWRITE,nid_fi)
	286	IF (ierr.NE.NF_NOERR) THEN
	287	write(6,*)' Problem opening file:',fichnom
	288	write(6,*)' ierr = ', ierr
	289	CALL ABORT
	290	ENDIF
	291
	292	tab0 = 0
	293	Lmodif = 0
	294
	295	endif
	296
	297	c-----------------------------------------------------------------------
	298	c Lecture du tableau des parametres du run (pour la dynamique)
	299	c-----------------------------------------------------------------------
	300
	301	if (choix_1.eq.0) then
	302
	303	write(,) 'reading tab_cntrl START_ARCHIVE'
	304	c
	305	ierr = NF_INQ_VARID (nid, "controle", nvarid)
	306	#ifdef NC_DOUBLE
	307	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tab_cntrl)
	308	#else
	309	ierr = NF_GET_VAR_REAL(nid, nvarid, tab_cntrl)
	310	#endif
	311	c
	312	else if (choix_1.eq.1) then
	313
	314	write(,) 'reading tab_cntrl START'
	315	c
	316	ierr = NF_INQ_VARID (nid_dyn, "controle", nvarid)
	317	#ifdef NC_DOUBLE
	318	ierr = NF_GET_VAR_DOUBLE(nid_dyn, nvarid, tab_cntrl)
	319	#else
	320	ierr = NF_GET_VAR_REAL(nid_dyn, nvarid, tab_cntrl)
	321	#endif
	322	c
	323	write(,) 'reading tab_cntrl STARTFI'
	324	c
	325	ierr = NF_INQ_VARID (nid_fi, "controle", nvarid)
	326	#ifdef NC_DOUBLE
	327	ierr = NF_GET_VAR_DOUBLE(nid_fi, nvarid, tab_cntrl_bis)
	328	#else
	329	ierr = NF_GET_VAR_REAL(nid_fi, nvarid, tab_cntrl_bis)
	330	#endif
	331	c
	332	do i=1,50
	333	tab_cntrl(i+50)=tab_cntrl_bis(i)
	334	enddo
	335	write(,) 'printing tab_cntrl', tab_cntrl
	336	do i=1,100
	337	write(,) i,tab_cntrl(i)
	338	enddo
	339
	340	endif
	341	c-----------------------------------------------------------------------
	342	c Initialisation des constantes dynamique
	343	c-----------------------------------------------------------------------
	344
	345	kappa = tab_cntrl(9)
	346	etot0 = tab_cntrl(12)
	347	ptot0 = tab_cntrl(13)
	348	ztot0 = tab_cntrl(14)
	349	stot0 = tab_cntrl(15)
	350	ang0 = tab_cntrl(16)
	351	write(,) "Newstart: kappa,etot0,ptot0,ztot0,stot0,ang0"
	352	write(,) kappa,etot0,ptot0,ztot0,stot0,ang0
	353
	354	c-----------------------------------------------------------------------
	355	c Lecture du tab_cntrl et initialisation des constantes physiques
	356	c - pour start: Lmodif = 0 => pas de modifications possibles
	357	c (modif dans le tabfi de readfi + loin)
	358	c - pour start_archive: Lmodif = 1 => modifications possibles
	359	c-----------------------------------------------------------------------
	360	if (choix_1.eq.0) then
[1130]	361	! tabfi requires that input file be first opened by open_startphy(fichnom)
	362	fichnom = 'start_archive.nc'
	363	call open_startphy(fichnom)
[38]	364	call tabfi (nid,Lmodif,tab0,day_ini,lllm,p_rad,
	365	. p_omeg,p_g,p_mugaz,p_daysec,time)
	366	else if (choix_1.eq.1) then
[1130]	367	fichnom = 'startfi.nc'
	368	call open_startphy(fichnom)
[38]	369	call tabfi (nid_fi,Lmodif,tab0,day_ini,lllm,p_rad,
	370	. p_omeg,p_g,p_mugaz,p_daysec,time)
	371	endif
	372
	373	rad = p_rad
	374	omeg = p_omeg
	375	g = p_g
	376	mugaz = p_mugaz
	377	daysec = p_daysec
	378
	379
	380	c=======================================================================
	381	c INITIALISATIONS DIVERSES
	382	c=======================================================================
	383
	384	day_step=180 !?! Note: day_step is a common in "control.h"
	385	CALL defrun_new( 99, .TRUE. )
[2825]	386	dtvr = daysec/REAL(day_step)
[38]	387	CALL iniconst
	388	CALL inigeom
	389	idum=-1
	390	idum=0
	391
[1543]	392	! Initialize the physics
	393	CALL iniphysiq(iim,jjm,llm,
	394	& (jjm-1)*iim+2,comm_lmdz,
	395	& daysec,day_ini,dtphys,
	396	& rlatu,rlatv,rlonu,rlonv,
	397	& aire,cu,cv,rad,g,r,cpp,
	398	& 1)
[38]	399
	400	c=======================================================================
	401	c lecture topographie, albedo, inertie thermique, relief sous-maille
	402	c=======================================================================
	403
	404	if (choix_1.ne.1) then ! pour ne pas avoir besoin du fichier
	405	! surface.dat dans le cas des start
	406
	407	c do while((relief(1:3).ne.'mol').AND.(relief(1:3).ne.'pla'))
	408	c write(,)
	409	c write(,) 'choix du relief (mola,pla)'
	410	c write(,) '(Topographie MGS MOLA, plat)'
	411	c read(*,fmt='(a3)') relief
	412	relief="mola"
	413	c enddo
	414
[224]	415	CALL datareadnc(relief,phis,alb,surfith,z0S,
[1974]	416	& zmeaS,zstdS,zsigS,zgamS,ztheS,
[2079]	417	& hmonsS,summitS,baseS,zavgS)
[38]	418
	419	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
	420	! CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ith,ithfi)
	421	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,surfith,surfithfi)
	422	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
[224]	423	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,z0S,z0)
[38]	424	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zmeaS,zmea)
	425	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zstdS,zstd)
	426	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zsigS,zsig)
	427	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zgamS,zgam)
	428	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ztheS,zthe)
[1974]	429	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,hmonsS,hmons)
[2079]	430	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,summitS,summit)
	431	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,baseS,base)
[1974]	432	! CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zavgS,zavg)
[38]	433
	434	endif ! of if (choix_1.ne.1)
	435
	436
	437	c=======================================================================
	438	c Lecture des fichiers (start ou start_archive)
	439	c=======================================================================
	440
	441	if (choix_1.eq.0) then
	442
	443	write(,) 'Reading file START_ARCHIVE'
[1047]	444	CALL lect_start_archive(ngridmx,llm,nqtot,
[2943]	445	& date,tsurf,tsoil,inertiesoil,albedo,emis,q2,
[2828]	446	& t,ucov,vcov,ps,teta,phisold_newgrid,q,qsurf,
[3025]	447	& tauscaling,totcloudfrac,surfith,nid,watercap,perenial_co2ice)
[38]	448	write(,) "OK, read start_archive file"
	449	! copy soil thermal inertia
	450	ithfi(:,:)=inertiedat(:,:)
	451
	452	ierr= NF_CLOSE(nid)
	453
	454	else if (choix_1.eq.1) then ! c'est l'appel a tabfi de phyeta0 qui
	455	! permet de changer les valeurs du
	456	! tab_cntrl Lmodif=1
	457	tab0=0
	458	Lmodif=1 ! Lmodif set to 1 to allow modifications in phyeta0
	459	write(,) 'Reading file START'
	460	fichnom = 'start.nc'
[1421]	461	CALL dynetat0(fichnom,vcov,ucov,teta,q,masse,
[1944]	462	& ps,phis,time)
[38]	463
	464	write(,) 'Reading file STARTFI'
	465	fichnom = 'startfi.nc'
[1047]	466	CALL phyetat0 (fichnom,tab0,Lmodif,nsoilmx,ngridmx,llm,nqtot,
[1944]	467	& day_ini,time,tsurf,tsoil,albedo,emis,
[2828]	468	& q2,qsurf,tauscaling,totcloudfrac,
[2999]	469	& wstar,watercap,perenial_co2ice,
	470	& def_slope,def_slope_mean,subslope_dist)
[38]	471
	472	! copy albedo and soil thermal inertia
	473	do i=1,ngridmx
	474	albfi(i) = albedodat(i)
	475	do j=1,nsoilmx
	476	ithfi(i,j) = inertiedat(i,j)
	477	enddo
	478	! build a surfithfi(:) using 1st layer of ithfi(:), which might
	479	! be neede later on if reinitializing soil thermal inertia
	480	surfithfi(i)=ithfi(i,1)
	481	enddo
	482
	483	else
	484	CALL exit(1)
	485	endif
	486
	487	dtvr = daysec/REAL(day_step)
	488	dtphys = dtvr * REAL(iphysiq)
	489
	490	c=======================================================================
	491	c
	492	c=======================================================================
	493	! If tracer names follow 'old' convention (q01, q02, ...) then
	494	! rename them
	495	count=0
[1036]	496	do iq=1,nqtot
[38]	497	txt=" "
	498	write(txt,'(a1,i2.2)') 'q',iq
[1130]	499	if (txt.eq.tname(iq)) then
[38]	500	count=count+1
	501	endif
[1036]	502	enddo ! of do iq=1,nqtot
[38]	503
[618]	504	! initialize tracer names noms(:) and indexes (igcm_co2, igcm_h2o_vap, ...)
[1224]	505	call initracer(ngridmx,nqtot,qsurf)
[618]	506
[1036]	507	if (count.eq.nqtot) then
[38]	508	write(,) 'Newstart: updating tracer names'
[1130]	509	! copy noms(:) to tname(:) to have matching tracer names in physics
[618]	510	! and dynamics
[1130]	511	tname(1:nqtot)=noms(1:nqtot)
[38]	512	endif
	513
	514	c=======================================================================
	515	c
	516	c=======================================================================
	517
	518	do ! infinite loop on list of changes
	519
	520	write(,)
	521	write(,)
	522	write(,) 'List of possible changes :'
[618]	523	write(,) '~~~~~~~~~~~~~~~~~~~~~~~~~~'
[38]	524	write(,)
[618]	525	write(,) 'flat : no topography ("aquaplanet")'
	526	write(,) 'bilball : uniform albedo and thermal inertia'
	527	write(,) 'z0 : set a uniform surface roughness length'
	528	write(,) 'coldspole : cold subsurface and high albedo at
	529	$ S.Pole'
	530	write(,) 'qname : change tracer name'
	531	write(,) 'q=0 : ALL tracer =zero'
[2217]	532	write(,) 'q=factor : change tracer value by a multiplicative
	533	& factor'
[618]	534	write(,) 'q=x : give a specific uniform value to one
	535	$ tracer'
	536	write(,) 'q=profile : specify a profile for a tracer'
[1130]	537	write(,) 'freedust : rescale dust to a true value'
[618]	538	write(,) 'ini_q : tracers initialization for chemistry
	539	$ and water vapour'
	540	write(,) 'ini_q-h2o : tracers initialization for chemistry
	541	$ only'
[1232]	542	write(,) 'composition : change atm main composition: CO2,N2,Ar,
	543	$ O2,CO'
[2312]	544	write(,) 'inihdo : initialize HDO'
[618]	545	write(,) 'ini_h2osurf : reinitialize surface water ice '
	546	write(,) 'noglacier : Remove tropical H2O ice if \|lat\|<45'
	547	write(,) 'watercapn : H20 ice on permanent N polar cap '
	548	write(,) 'watercaps : H20 ice on permanent S polar cap '
	549	write(,) 'wetstart : start with a wet atmosphere'
	550	write(,) 'isotherm : Isothermal Temperatures, wind set to
	551	$ zero'
[2828]	552	write(,) 'co2ice=0 : remove CO2 polar cap i.e.
	553	$ qsurf(co2)=0 '
[618]	554	write(,) 'ptot : change total pressure'
	555	write(,) 'therm_ini_s : set soil thermal inertia to reference
	556	$ surface values'
	557	write(,) 'subsoilice_n : put deep underground ice layer in
	558	$ northern hemisphere'
	559	write(,) 'subsoilice_s : put deep underground ice layer in
	560	$ southern hemisphere'
	561	write(,) 'mons_ice : put underground ice layer according
	562	$ to MONS derived data'
[2909]	563	write(,) 'nslope : Change the number of subgrid scale
	564	$ slope'
[38]	565
	566	write(,)
	567	write(,) 'Change to perform ?'
	568	write(,) ' (enter keyword or return to end)'
	569	write(,)
	570
	571	read(*,fmt='(a20)') modif
	572	if (modif(1:1) .eq. ' ') exit ! exit loop on changes
	573
	574	write(,)
[164]	575	write(,) trim(modif) , ' : '
[38]	576
	577	c 'flat : no topography ("aquaplanet")'
	578	c -------------------------------------
[164]	579	if (trim(modif) .eq. 'flat') then
[38]	580	c set topo to zero
[1403]	581	phis(:,:)=0
[38]	582	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
	583	write(,) 'topography set to zero.'
	584	write(,) 'WARNING : the subgrid topography parameters',
	585	& ' were not set to zero ! => set calllott to F'
	586
	587	c Choice for surface pressure
	588	yes=' '
	589	do while ((yes.ne.'y').and.(yes.ne.'n'))
	590	write(,) 'Do you wish to choose homogeneous surface',
	591	& 'pressure (y) or let newstart interpolate ',
	592	& ' the previous field (n)?'
	593	read(*,fmt='(a)') yes
	594	end do
	595	if (yes.eq.'y') then
	596	flagps0=.true.
	597	write(,) 'New value for ps (Pa) ?'
	598	201 read(,,iostat=ierr) patm
	599	if(ierr.ne.0) goto 201
	600	write(,)
	601	write(,) ' new ps everywhere (Pa) = ', patm
	602	write(,)
	603	do j=1,jjp1
	604	do i=1,iip1
	605	ps(i,j)=patm
	606	enddo
	607	enddo
	608	end if
	609
	610	c bilball : albedo, inertie thermique uniforme
	611	c --------------------------------------------
[164]	612	else if (trim(modif) .eq. 'bilball') then
[38]	613	write(,) 'constante albedo and iner.therm:'
	614	write(,) 'New value for albedo (ex: 0.25) ?'
	615	101 read(,,iostat=ierr) alb_bb
	616	if(ierr.ne.0) goto 101
	617	write(,)
	618	write(,) ' uniform albedo (new value):',alb_bb
	619	write(,)
	620
	621	write(,) 'New value for thermal inertia (eg: 247) ?'
	622	102 read(,,iostat=ierr) ith_bb
	623	if(ierr.ne.0) goto 102
	624	write(,) 'uniform thermal inertia (new value):',ith_bb
	625	DO j=1,jjp1
	626	DO i=1,iip1
	627	alb(i,j) = alb_bb ! albedo
	628	do isoil=1,nsoilmx
	629	ith(i,j,isoil) = ith_bb ! thermal inertia
	630	enddo
	631	END DO
	632	END DO
	633	! CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ith,ithfi)
	634	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
	635	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
[224]	636
	637	! also reset surface roughness length to default value
	638	write(,) 'surface roughness length set to:',z0_default,' m'
	639	z0(:)=z0_default
[38]	640
[224]	641	! z0 : set surface roughness length to a constant value
	642	! -----------------------------------------------------
	643	else if (trim(modif) .eq. 'z0') then
	644	write(,) 'set a uniform surface roughness length'
	645	write(,) ' value for z0_default (ex: ',z0_default,')?'
	646	ierr=1
	647	do while (ierr.ne.0)
	648	read(,,iostat=ierr) z0_default
	649	enddo
	650	z0(:)=z0_default
	651
[38]	652	c coldspole : sous-sol de la calotte sud toujours froid
	653	c -----------------------------------------------------
[164]	654	else if (trim(modif) .eq. 'coldspole') then
[38]	655	write(,)'new value for the subsurface temperature',
	656	& ' beneath the permanent southern polar cap ? (eg: 141 K)'
	657	103 read(,,iostat=ierr) tsud
	658	if(ierr.ne.0) goto 103
	659	write(,)
	660	write(,) ' new value of the subsurface temperature:',tsud
	661	c nouvelle temperature sous la calotte permanente
[2909]	662	do islope=1,nslope
	663	do l=2,nsoilmx
	664	tsoil(ngridmx,l,islope) = tsud
	665	end do
	666	enddo
[38]	667
	668
	669	write(,)'new value for the albedo',
	670	& 'of the permanent southern polar cap ? (eg: 0.75)'
	671	104 read(,,iostat=ierr) albsud
	672	if(ierr.ne.0) goto 104
	673	write(,)
	674
	675	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	676	c Option 1: only the albedo of the pole is modified :
	677	albfi(ngridmx)=albsud
	678	write(,) 'ig=',ngridmx,' albedo perennial cap ',
	679	& albfi(ngridmx)
	680
	681	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	682	c Option 2 A haute resolution : coordonnee de la vrai calotte ~
	683	c DO j=1,jjp1
	684	c DO i=1,iip1
	685	c ig=1+(j-2)*iim +i
	686	c if(j.eq.1) ig=1
	687	c if(j.eq.jjp1) ig=ngridmx
	688	c if ((rlatu(j)*180./pi.lt.-84.).and.
	689	c & (rlatu(j)*180./pi.gt.-91.).and.
	690	c & (rlonv(i)*180./pi.gt.-91.).and.
	691	c & (rlonv(i)*180./pi.lt.0.)) then
	692	cc albedo de la calotte permanente fixe a albsud
	693	c alb(i,j)=albsud
	694	c write(,) 'lat=',rlatu(j)*180./pi,
	695	c & ' lon=',rlonv(i)*180./pi
	696	cc fin de la condition sur les limites de la calotte permanente
	697	c end if
	698	c ENDDO
	699	c ENDDO
	700	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	701
	702	c CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
	703
	704
	705	c ptot : Modification of the total pressure: ice + current atmosphere
	706	c -------------------------------------------------------------------
[164]	707	else if (trim(modif) .eq. 'ptot') then
[38]	708
	709	c calcul de la pression totale glace + atm actuelle
	710	patm=0.
	711	airetot=0.
	712	pcap=0.
	713	DO j=1,jjp1
	714	DO i=1,iim
	715	ig=1+(j-2)*iim +i
	716	if(j.eq.1) ig=1
	717	if(j.eq.jjp1) ig=ngridmx
	718	patm = patm + ps(i,j)*aire(i,j)
	719	airetot= airetot + aire(i,j)
[2909]	720	DO islope=1,nslope
	721	pcap = pcap + aire(i,j)qsurf(ig,igcm_co2,islope)g
	722	& subslope_dist(ig,islope)/cos(def_slope_mean(islope)pi/180.)
	723	ENDDO
[38]	724	ENDDO
	725	ENDDO
	726	ptoto = pcap + patm
	727
	728	print*,'Current total pressure at surface (co2 ice + atm) ',
	729	& ptoto/airetot
	730
	731	print*,'new value?'
	732	read(,) ptotn
	733	ptotn=ptotn*airetot
	734	patmn=ptotn-pcap
	735	print*,'ptoto,patm,ptotn,patmn'
	736	print*,ptoto,patm,ptotn,patmn
	737	print*,'Mult. factor for pressure (atm only)', patmn/patm
	738	do j=1,jjp1
	739	do i=1,iip1
	740	ps(i,j)=ps(i,j)*patmn/patm
	741	enddo
	742	enddo
	743
	744	c Correction pour la conservation des traceurs
	745	yes=' '
	746	do while ((yes.ne.'y').and.(yes.ne.'n'))
	747	write(,) 'Do you wish to conserve tracer total mass (y)',
	748	& ' or tracer mixing ratio (n) ?'
	749	read(*,fmt='(a)') yes
	750	end do
	751
	752	if (yes.eq.'y') then
	753	write(,) 'OK : conservation of tracer total mass'
[1036]	754	DO iq =1, nqtot
[38]	755	DO l=1,llm
	756	DO j=1,jjp1
	757	DO i=1,iip1
	758	q(i,j,l,iq)=q(i,j,l,iq)*patm/patmn
	759	ENDDO
	760	ENDDO
	761	ENDDO
	762	ENDDO
	763	else
	764	write(,) 'OK : conservation of tracer mixing ratio'
	765	end if
	766
	767	c qname : change tracer name
	768	c --------------------------
	769	else if (trim(modif).eq.'qname') then
	770	yes='y'
	771	do while (yes.eq.'y')
	772	write(,) 'Which tracer name do you want to change ?'
[1036]	773	do iq=1,nqtot
[1130]	774	write(*,'(i3,a3,a20)')iq,' : ',trim(tname(iq))
[38]	775	enddo
	776	write(*,'(a35,i3)')
[1036]	777	& '(enter tracer number; between 1 and ',nqtot
[38]	778	write(,)' or any other value to quit this option)'
	779	read(,) iq
[1036]	780	if ((iq.ge.1).and.(iq.le.nqtot)) then
[1130]	781	write(,)'Change tracer name ',trim(tname(iq)),' to ?'
[38]	782	read(,) txt
[1130]	783	tname(iq)=txt
[38]	784	write(,)'Do you want to change another tracer name (y/n)?'
	785	read(*,'(a)') yes
	786	else
	787	! inapropiate value of iq; quit this option
	788	yes='n'
[1036]	789	endif ! of if ((iq.ge.1).and.(iq.le.nqtot))
[38]	790	enddo ! of do while (yes.ne.'y')
	791
	792	c q=0 : set tracers to zero
	793	c -------------------------
[164]	794	else if (trim(modif) .eq. 'q=0') then
[38]	795	c mise a 0 des q (traceurs)
	796	write(,) 'Tracers set to 0 (1.E-30 in fact)'
[1036]	797	DO iq =1, nqtot
[38]	798	DO l=1,llm
	799	DO j=1,jjp1
	800	DO i=1,iip1
	801	q(i,j,l,iq)=1.e-30
	802	ENDDO
	803	ENDDO
	804	ENDDO
	805	ENDDO
	806
	807	c set surface tracers to zero
[1036]	808	DO iq =1, nqtot
[38]	809	DO ig=1,ngridmx
[2909]	810	DO islope=1,nslope
	811	qsurf(ig,iq,islope)=0.
	812	ENDDO
[38]	813	ENDDO
	814	ENDDO
	815
[2217]	816	c q=factor : change value of tracer by a multiplicative factor
	817	c ------------------------------------------------------------
	818	else if (trim(modif) .eq. 'q=factor') then
	819	write(,) 'Which tracer do you want to modify ?'
	820	do iq=1,nqtot
	821	write(,)iq,' : ',trim(tname(iq))
	822	enddo
	823	write(,) '(choose between 1 and ',nqtot,')'
	824	read(,) iq
	825	if ((iq.lt.1).or.(iq.gt.nqtot)) then
	826	! wrong value for iq, go back to menu
	827	write(,) "wrong input value:",iq
	828	cycle
	829	endif
	830	write(,)"factor to multiply current mixing ratio by?"
	831	read(,) val
	832
	833	q(1:iip1,1:jjp1,1:llm,iq)=q(1:iip1,1:jjp1,1:llm,iq)*val
[2909]	834	qsurf(1:ngridmx,iq,:)=qsurf(1:ngridmx,iq,:)*val
[2217]	835
[38]	836	c q=x : initialise tracer manually
	837	c --------------------------------
[164]	838	else if (trim(modif) .eq. 'q=x') then
[38]	839	write(,) 'Which tracer do you want to modify ?'
[1036]	840	do iq=1,nqtot
[1130]	841	write(,)iq,' : ',trim(tname(iq))
[38]	842	enddo
[1036]	843	write(,) '(choose between 1 and ',nqtot,')'
[38]	844	read(,) iq
[1036]	845	if ((iq.lt.1).or.(iq.gt.nqtot)) then
[171]	846	! wrong value for iq, go back to menu
	847	write(,) "wrong input value:",iq
	848	cycle
	849	endif
[1130]	850	write(,)'mixing ratio of tracer ',trim(tname(iq)),
[38]	851	& ' ? (kg/kg)'
	852	read(,) val
	853	DO l=1,llm
	854	DO j=1,jjp1
	855	DO i=1,iip1
	856	q(i,j,l,iq)=val
	857	ENDDO
	858	ENDDO
	859	ENDDO
[1130]	860	write(,) 'SURFACE value of tracer ',trim(tname(iq)),
[38]	861	& ' ? (kg/m2)'
	862	read(,) val
	863	DO ig=1,ngridmx
[2909]	864	DO islope=1,nslope
	865	qsurf(ig,iq,islope)=val
	866	ENDDO
[38]	867	ENDDO
	868
[563]	869	c q=profile : initialize tracer with a given profile
	870	c --------------------------------------------------
	871	else if (trim(modif) .eq. 'q=profile') then
	872	write(,) 'Tracer profile will be sought in ASCII file'
	873	write(,) "'profile_tracer' where 'tracer' is tracer name"
	874	write(,) "(one value per line in file; starting with"
	875	write(,) "surface value, the 1st atmospheric layer"
	876	write(,) "followed by 2nd, etc. up to top of atmosphere)"
	877	write(,) 'Which tracer do you want to set?'
[1036]	878	do iq=1,nqtot
[1130]	879	write(,)iq,' : ',trim(tname(iq))
[563]	880	enddo
[1036]	881	write(,) '(choose between 1 and ',nqtot,')'
[563]	882	read(,) iq
[1036]	883	if ((iq.lt.1).or.(iq.gt.nqtot)) then
[563]	884	! wrong value for iq, go back to menu
	885	write(,) "wrong input value:",iq
	886	cycle
	887	endif
	888	! look for input file 'profile_tracer'
[1130]	889	txt="profile_"//trim(tname(iq))
[563]	890	open(41,file=trim(txt),status='old',form='formatted',
	891	& iostat=ierr)
	892	if (ierr.eq.0) then
	893	! OK, found file 'profile_...', load the profile
	894	do l=1,llm+1
	895	read(41,*,iostat=ierr) profile(l)
	896	if (ierr.ne.0) then ! something went wrong
	897	exit ! quit loop
	898	endif
	899	enddo
	900	if (ierr.eq.0) then
	901	! initialize tracer values
[2909]	902	do islope=1,nslope
	903	qsurf(:,iq,islope)=profile(1)
	904	enddo
[563]	905	do l=1,llm
	906	q(:,:,l,iq)=profile(l+1)
	907	enddo
[1130]	908	write(,)'OK, tracer ',trim(tname(iq)),
	909	& ' initialized ','using values from file ',trim(txt)
[563]	910	else
	911	write(,)'problem reading file ',trim(txt),' !'
[1130]	912	write(,)'No modifications to tracer ',trim(tname(iq))
[563]	913	endif
	914	else
	915	write(,)'Could not find file ',trim(txt),' !'
[1130]	916	write(,)'No modifications to tracer ',trim(tname(iq))
[563]	917	endif
	918
[1208]	919	c convert dust from virtual to true values
[1130]	920	c --------------------------------------------------
	921	else if (trim(modif) .eq. 'freedust') then
[1208]	922	if (minval(tauscaling) .lt. 0) then
	923	write(,) 'WARNING conversion factor negative'
	924	write(,) 'This is probably because it was not present
	925	&in the file'
	926	write(,) 'A constant conversion is used instead.'
	927	tauscaling(:) = 1.e-3
	928	endif
	929	CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,tauscaling,tauscadyn)
[1130]	930	do l=1,llm
	931	do j=1,jjp1
	932	do i=1,iip1
	933	if (igcm_dust_number .ne. 0)
[1208]	934	& q(i,j,l,igcm_dust_number) =
	935	& q(i,j,l,igcm_dust_number) * tauscadyn(i,j)
[1130]	936	if (igcm_dust_mass .ne. 0)
[1208]	937	& q(i,j,l,igcm_dust_mass) =
	938	& q(i,j,l,igcm_dust_mass) * tauscadyn(i,j)
[1130]	939	if (igcm_ccn_number .ne. 0)
[1208]	940	& q(i,j,l,igcm_ccn_number) =
	941	& q(i,j,l,igcm_ccn_number) * tauscadyn(i,j)
[1130]	942	if (igcm_ccn_mass .ne. 0)
[1208]	943	& q(i,j,l,igcm_ccn_mass) =
	944	& q(i,j,l,igcm_ccn_mass) * tauscadyn(i,j)
[1130]	945	end do
	946	end do
	947	end do
[563]	948
[1208]	949	tauscaling(:) = 1.
	950
[1130]	951	! We want to have the very same value at lon -180 and lon 180
	952	do l = 1,llm
	953	do j = 1,jjp1
	954	do iq = 1,nqtot
	955	q(iip1,j,l,iq) = q(1,j,l,iq)
	956	end do
	957	end do
	958	end do
	959
[1208]	960	write(,) 'done rescaling to true vale'
	961
[38]	962	c ini_q : Initialize tracers for chemistry
	963	c -----------------------------------------------
[164]	964	else if (trim(modif) .eq. 'ini_q') then
[618]	965	flagh2o = 1
	966	flagthermo = 0
	967	yes=' '
[38]	968	c For more than 32 layers, possible to initiate thermosphere only
	969	if (llm.gt.32) then
	970	do while ((yes.ne.'y').and.(yes.ne.'n'))
	971	write(,)'',
	972	& 'initialisation for thermosphere only? (y/n)'
	973	read(*,fmt='(a)') yes
	974	if (yes.eq.'y') then
[618]	975	flagthermo=1
[38]	976	else
[618]	977	flagthermo=0
[38]	978	endif
	979	enddo
	980	endif
	981
[1231]	982	call inichim_newstart(ngridmx, nqtot, q, qsurf, ps,
[1047]	983	& flagh2o, flagthermo)
[664]	984
	985	! We want to have the very same value at lon -180 and lon 180
	986	do l = 1,llm
	987	do j = 1,jjp1
[1036]	988	do iq = 1,nqtot
[664]	989	q(iip1,j,l,iq) = q(1,j,l,iq)
	990	end do
	991	end do
	992	end do
	993
[618]	994	write(,) 'inichim_newstart: chemical species and
	995	$ water vapour initialised'
[38]	996
[664]	997	c ini_q-h2o : as above except for the water vapour tracer
[38]	998	c ------------------------------------------------------
[618]	999	else if (trim(modif) .eq. 'ini_q-h2o') then
	1000	flagh2o = 0
	1001	flagthermo = 0
	1002	yes=' '
[38]	1003	! for more than 32 layers, possible to initiate thermosphere only
	1004	if(llm.gt.32) then
	1005	do while ((yes.ne.'y').and.(yes.ne.'n'))
	1006	write(,)'',
	1007	& 'initialisation for thermosphere only? (y/n)'
	1008	read(*,fmt='(a)') yes
	1009	if (yes.eq.'y') then
[618]	1010	flagthermo=1
[38]	1011	else
[618]	1012	flagthermo=0
[38]	1013	endif
	1014	enddo
	1015	endif
[664]	1016
[1231]	1017	call inichim_newstart(ngridmx, nqtot, q, qsurf, ps,
	1018	& flagh2o, flagthermo)
[664]	1019
	1020	! We want to have the very same value at lon -180 and lon 180
	1021	do l = 1,llm
	1022	do j = 1,jjp1
[1036]	1023	do iq = 1,nqtot
[664]	1024	q(iip1,j,l,iq) = q(1,j,l,iq)
	1025	end do
	1026	end do
	1027	end do
	1028
[618]	1029	write(,) 'inichim_newstart: chemical species initialised
	1030	$ (except water vapour)'
[38]	1031
[2312]	1032	c inihdo : initialize HDO with user D/H value
	1033	c --------------------------------------------------------
	1034	else if (trim(modif) .eq. 'inihdo') then
	1035	! check that there is indeed a water vapor tracer
	1036	if (igcm_h2o_vap.eq.0) then
	1037	write(,) "No water vapour tracer! Can't use this option"
	1038	stop
	1039	endif
	1040
	1041	write(,)'Input D/H ratio (in SMOW)'
	1042	write(,)'If value is <0 then HDO=H2O'
	1043	303 read(,, iostat=ierr) DoverH
	1044	if(ierr.ne.0) goto 303
	1045
	1046	DoverH = DoverH * 2 * 155.76e-6
	1047	if (DoverH.lt.0.0) then
	1048	DoverH = 1.
	1049	endif
	1050	!D/H (SMOW) = 155.76e-6 so HDO/H2O is twice that
	1051
	1052	do ig=1,ngridmx
[2909]	1053	do islope=1,nslope
	1054	qsurf(ig,igcm_h2o_ice,islope)=
	1055	& max(0.,qsurf(ig,igcm_h2o_ice,islope))
	1056	enddo
[2312]	1057	end do
	1058
[2378]	1059	! Update the hdo tracers
	1060	q(1:iip1,1:jjp1,1:llm,igcm_hdo_vap)
	1061	& =q(1:iip1,1:jjp1,1:llm,igcm_h2o_vap)* DoverH
	1062	q(1:iip1,1:jjp1,1:llm,igcm_hdo_ice)
	1063	& =q(1:iip1,1:jjp1,1:llm,igcm_h2o_ice)* DoverH
[2312]	1064
[2909]	1065	qsurf(1:ngridmx,igcm_hdo_ice,:)
	1066	& =qsurf(1:ngridmx,igcm_h2o_ice,:)*DoverH
[2312]	1067
[2378]	1068
	1069
[1232]	1070	c composition : change main composition: CO2,N2,Ar,O2,CO (FF 03/2014)
	1071	c --------------------------------------------------------
	1072	else if (trim(modif) .eq. 'composition') then
	1073	write(,) "Lat (degN) lon (degE) of the reference site ?"
	1074	write(,) "e.g. MSL : -4.5 137. "
	1075	301 read(,,iostat=ierr) latref, lonref
	1076	if(ierr.ne.0) goto 301
	1077
	1078
	1079	! Select GCM point close to reference site
	1080	dlonmin =90.
	1081	DO i=1,iip1-1
	1082	if (abs(rlonv(i)*180./pi -lonref).lt.dlonmin)then
	1083	iref=i
	1084	dlonmin=abs(rlonv(i)*180./pi -lonref)
	1085	end if
	1086	ENDDO
	1087	dlatmin =45.
	1088	DO j=1,jjp1
	1089	if (abs(rlatu(j)*180./pi -latref).lt.dlatmin)then
	1090	jref=j
	1091	dlatmin=abs(rlatu(j)*180./pi -latref)
	1092	end if
	1093	ENDDO
	1094	write(,) "In GCM : lat= " , rlatu(jref)*180./pi
	1095	write(,) "In GCM : lon= " , rlonv(iref)*180./pi
	1096	write(,)
	1097
	1098	! Compute air molar mass at reference site
[2175]	1099	Smmr=0.
	1100	Sn = 0.
	1101	write(,) 'igcm_co2 = ', igcm_co2
	1102	write(,) 'igcm_n2 = ', igcm_n2
	1103	write(,) 'igcm_ar = ', igcm_ar
	1104	write(,) 'igcm_o2 = ', igcm_o2
	1105	write(,) 'igcm_co = ', igcm_co
	1106	write(,) noms
[1232]	1107	do iq=1,nqtot
	1108	if ((iq.eq.igcm_co2).or.(iq.eq.igcm_n2)
	1109	& .or. (iq.eq.igcm_ar).or.(iq.eq.igcm_o2)
	1110	& .or. (iq.eq.igcm_co)) then
	1111	Smmr=Smmr+q(iref,jref,1,iq)
	1112	Sn=Sn+q(iref,jref,1,iq)/mmol(iq)
	1113	end if
	1114	end do
[2175]	1115	! Special case : implicit non-co2 gases ! JN 11/2019
	1116	if ((igcm_n2.eq.0) .or. (igcm_ar.eq.0)) then
	1117	write(,) "Warning : non-co2 gases are implicit : "
	1118	write(,) "At reference site : "
	1119	! write(,) "q= ", q(iref, jref, 1,igcm_co2)
	1120	write(,) "Sum of mass mix. ratio (ie MMR(co2))=",Smmr
	1121	Mair_old = 44.0Smmr + 33.87226017157708(1-Smmr)
	1122
	1123	! 33.87226017157708 is the
	1124	! molar mass of non-co2 atmosphere measured by MSL at Ls ~184
	1125
[2217]	1126	else
	1127	! Assume co2/n2/ar/o2/co are available
	1128	Mair_old=(q(iref,jref,1,igcm_co2)*mmol(igcm_co2)
	1129	& +q(iref,jref,1,igcm_n2)*mmol(igcm_n2)
	1130	& +q(iref,jref,1,igcm_ar)*mmol(igcm_ar)
	1131	& +q(iref,jref,1,igcm_o2)*mmol(igcm_o2)
	1132	& +q(iref,jref,1,igcm_co)*mmol(igcm_co))/Smmr
[2175]	1133	end if
[1232]	1134
[2217]	1135	write(,)
	1136	& "Air molar mass (g/mol) at reference site= ",Mair_old
	1137
[1232]	1138	! Ask for new volume mixing ratio at reference site
	1139	Svmr =0.
	1140	Sn =0.
[2175]	1141	coefvmr(igcm_co2)=1.
	1142
[1232]	1143	do iq=1,nqtot
	1144	coefvmr(iq) = 1.
	1145	if ((iq.eq.igcm_n2).or.(iq.eq.igcm_ar)
	1146	& .or. (iq.eq.igcm_o2).or.(iq.eq.igcm_co)) then
	1147
	1148	vmr_old=q(iref,jref,1,iq)*Mair_old/mmol(iq)
	1149	write(,) "Previous vmr("//trim(tname(iq))//")= ", vmr_old
	1150
	1151	if (iq.eq.igcm_n2) then
[2258]	1152	write(,) "New vmr(n2)? (MSL: 2.8e-02 at Ls~180,",
	1153	& " Trainer et al. 2019)"
[1232]	1154	endif
	1155	if (iq.eq.igcm_ar) then
[2258]	1156	write(,) "New vmr(ar)? (MSL: 2.1e-02 at Ls~180)"
[1232]	1157	endif
	1158	if (iq.eq.igcm_o2) then
[2258]	1159	write(,) "New vmr(o2)? (MSL: 1.7e-03 at Ls~180)"
[1232]	1160	endif
	1161	if (iq.eq.igcm_co) then
[2258]	1162	write(,) "New vmr(co)? (ACS: 1e-03 at Ls~180)"
[1232]	1163	endif
	1164	302 read(,,iostat=ierr) vmr_new
	1165	if(ierr.ne.0) goto 302
	1166	write(,) "New vmr("//trim(tname(iq))//")= ",vmr_new
	1167	write(,)
	1168	coefvmr(iq) = vmr_new/vmr_old
	1169	Svmr=Svmr+vmr_new
	1170	Sn=Sn+vmr_new*mmol(iq)
	1171	end if
	1172	enddo ! of do iq=1,nqtot
[2175]	1173
	1174	! Special case : implicit non-co2 gases JN 11/2019
	1175	if ((igcm_n2.eq.0) .or. (igcm_ar.eq.0)) then
	1176	write(,) "Warning : non-co2 gases are implicit"
	1177	vmr_old=q(iref,jref,1,igcm_co2)*Mair_old/mmol(igcm_co2)
	1178	write(,) "Previous vmr(co2)=", vmr_old
[2258]	1179	write(,) "New vmr(co2) ? (MSL: 0.947 at Ls~180)",
	1180	& " Trainer et al. 2019)"
[2175]	1181	666 read(,,iostat=ierr) vmr_new
	1182	if(ierr.ne.0) goto 666
	1183	coefvmr(igcm_co2) = vmr_new/vmr_old
	1184	Svmr=Svmr+vmr_new
	1185	Sn=vmr_new*mmol(igcm_co2) + (1-vmr_new)
	1186	& *33.87226017157708 ! Molar mass of non-co2 atm from MSL
	1187	end if
[1232]	1188	! Estimation of new Air molar mass at reference site (assuming vmr_co2 = 1-Svmr)
	1189	Mair_new = Sn + (1-Svmr)*mmol(igcm_co2)
[2175]	1190	! Estimation of new Air molar mass when non-co2 gases are implicit
	1191	if ((igcm_n2.eq.0) .or. (igcm_ar.eq.0)) then
	1192	Mair_new=vmr_new*mmol(igcm_co2) + (1-vmr_new)
	1193	& *33.87226017157708 ! Molar mass of non-co2 atm from MSL
	1194	write(,)
	1195	& "We consider non-co2 gases vmr measured from Curiosity"
	1196	end if
[1232]	1197	write(,)
	1198	& "NEW Air molar mass (g/mol) at reference site= ",Mair_new
	1199
	1200	! Compute mass mixing ratio changes
	1201	do iq=1,nqtot
[2175]	1202	if ((iq.eq.igcm_co2).or.(iq.eq.igcm_n2).or.(iq.eq.igcm_ar)
[1232]	1203	& .or. (iq.eq.igcm_o2).or.(iq.eq.igcm_co)) then
	1204	write(,) "Everywhere mmr("//trim(tname(iq))//
	1205	& ") is multiplied by ",coefvmr(iq)*Mair_old/Mair_new
	1206	end if
	1207	end do
	1208
[1390]	1209	! Recompute mass mixing ratios everywhere, and adjust mmr of most abundant species
	1210	! to keep sum of mmr constant.
[1232]	1211	do l=1,llm
	1212	do j=1,jjp1
	1213	do i=1,iip1
	1214	Smmr_old = 0.
	1215	Smmr_new = 0.
	1216	do iq=1,nqtot
[2175]	1217	if ((iq.eq.igcm_co2).or.(iq.eq.igcm_n2)
	1218	& .or.(iq.eq.igcm_ar)
[2258]	1219	& .or. (iq.eq.igcm_o2).or.(iq.eq.igcm_co)
	1220	& .or. (iq.eq.igcm_o) .or. (iq.eq. igcm_h2) ) then
[1232]	1221	Smmr_old = Smmr_old + q(i,j,l,iq) ! sum of old mmr
	1222	q(i,j,l,iq)=q(i,j,l,iq)coefvmr(iq)Mair_old/Mair_new
	1223	Smmr_new = Smmr_new + q(i,j,l,iq) ! sum of new mmr
	1224	end if
	1225	enddo
[2258]	1226	!iloc = maxloc(q(i,j,l,:))
	1227	iqmax=0 ; maxq=0
	1228	do iq=1,nqtot
	1229	if ((iq.eq.igcm_co2).or.(iq.eq.igcm_n2)
	1230	& .or.(iq.eq.igcm_ar)
	1231	& .or. (iq.eq.igcm_o2).or.(iq.eq.igcm_co)
	1232	& .or. (iq.eq.igcm_o) .or. (iq.eq. igcm_h2) ) then
	1233	if (q(i,j,l,iq).gt.maxq) then
	1234	maxq=q(i,j,l,iq)
	1235	iqmax=iq
	1236	endif
	1237	endif
	1238	enddo
	1239	!iqmax = iloc(1)
[1390]	1240	q(i,j,l,iqmax) = q(i,j,l,iqmax) + Smmr_old - Smmr_new
[1232]	1241	enddo
	1242	enddo
	1243	enddo
	1244
	1245	write(,)
[1390]	1246	& "The most abundant species is modified everywhere to keep "//
	1247	& "sum of mmr constant"
[1232]	1248	write(,) 'At reference site vmr(CO2)=',
	1249	& q(iref,jref,1,igcm_co2)*Mair_new/mmol(igcm_co2)
[2258]	1250	write(,) "Compared to MSL observation: vmr(CO2)= 0.947 "//
	1251	& "at Ls=180"
[1232]	1252
[2258]	1253	Sn = q(iref,jref,1,igcm_co2)*Mair_new/mmol(igcm_co2)
	1254	& + q(iref,jref,1,igcm_n2)*Mair_new/mmol(igcm_n2)
	1255	& + q(iref,jref,1,igcm_ar)*Mair_new/mmol(igcm_ar)
	1256	& + q(iref,jref,1,igcm_o2)*Mair_new/mmol(igcm_o2)
	1257	& + q(iref,jref,1,igcm_co)*Mair_new/mmol(igcm_co)
	1258
	1259	write(,) 'Sum of volume mixing ratios = ', Sn
	1260
[38]	1261	c wetstart : wet atmosphere with a north to south gradient
	1262	c --------------------------------------------------------
[164]	1263	else if (trim(modif) .eq. 'wetstart') then
[38]	1264	! check that there is indeed a water vapor tracer
	1265	if (igcm_h2o_vap.eq.0) then
	1266	write(,) "No water vapour tracer! Can't use this option"
	1267	stop
	1268	endif
	1269	DO l=1,llm
	1270	DO j=1,jjp1
[321]	1271	DO i=1,iip1-1
[38]	1272	q(i,j,l,igcm_h2o_vap)=150.e-6 * (rlatu(j)+pi/2.) / pi
	1273	ENDDO
[321]	1274	! We want to have the very same value at lon -180 and lon 180
	1275	q(iip1,j,l,igcm_h2o_vap) = q(1,j,l,igcm_h2o_vap)
[38]	1276	ENDDO
	1277	ENDDO
	1278
	1279	write(,) 'Water mass mixing ratio at north pole='
	1280	* ,q(1,1,1,igcm_h2o_vap)
	1281	write(,) '---------------------------south pole='
	1282	* ,q(1,jjp1,1,igcm_h2o_vap)
	1283
[480]	1284	c ini_h2osurf : reinitialize surface water ice
	1285	c --------------------------------------------------
	1286	else if (trim(modif) .eq. 'ini_h2osurf') then
	1287	write(,)'max surface ice left?(e.g. 0.2 kg/m2=200microns)'
	1288	207 read(,,iostat=ierr) val
	1289	if(ierr.ne.0) goto 207
	1290	write(,)'also set negative values of surf ice to 0'
	1291	do ig=1,ngridmx
[2909]	1292	do islope=1,nslope
	1293	qsurf(ig,igcm_h2o_ice,islope)=
	1294	& min(val,qsurf(ig,igcm_h2o_ice,islope))
	1295	qsurf(ig,igcm_h2o_ice,islope)=
	1296	& max(0.,qsurf(ig,igcm_h2o_ice,islope))
	1297	enddo
[480]	1298	end do
	1299
[38]	1300	c noglacier : remove tropical water ice (to initialize high res sim)
	1301	c --------------------------------------------------
[164]	1302	else if (trim(modif) .eq. 'noglacier') then
[38]	1303	do ig=1,ngridmx
	1304	j=(ig-2)/iim +2
	1305	if(ig.eq.1) j=1
	1306	write(,) 'OK: remove surface ice for \|lat\|<45'
	1307	if (abs(rlatu(j)*180./pi).lt.45.) then
[2909]	1308	do islope=1,nslope
	1309	qsurf(ig,igcm_h2o_ice,islope)=0.
	1310	enddo
[38]	1311	end if
	1312	end do
	1313
	1314
	1315	c watercapn : H20 ice on permanent northern cap
	1316	c --------------------------------------------------
[164]	1317	else if (trim(modif) .eq. 'watercapn') then
[38]	1318	do ig=1,ngridmx
	1319	j=(ig-2)/iim +2
	1320	if(ig.eq.1) j=1
	1321	if (rlatu(j)*180./pi.gt.80.) then
[2909]	1322	do islope=1,nslope
	1323	qsurf(ig,igcm_h2o_ice,islope)=1.e5
	1324	write(,) 'ig=',ig,', islope=', islope,
	1325	& ' H2O ice mass (kg/m2)= ',
	1326	& qsurf(ig,igcm_h2o_ice,islope)
[38]	1327	write(,)' ==> Ice mesh South boundary (deg)= ',
	1328	& rlatv(j)*180./pi
[2909]	1329	enddo
[38]	1330	end if
	1331	enddo
	1332
	1333	c watercaps : H20 ice on permanent southern cap
	1334	c -------------------------------------------------
[164]	1335	else if (trim(modif) .eq. 'watercaps') then
[38]	1336	do ig=1,ngridmx
	1337	j=(ig-2)/iim +2
	1338	if(ig.eq.1) j=1
	1339	if (rlatu(j)*180./pi.lt.-80.) then
[2909]	1340	do islope=1,nslope
	1341	qsurf(ig,igcm_h2o_ice,islope)=1.e5
	1342	write(,) 'ig=',ig,', islope=', islope,
	1343	& ' H2O ice mass (kg/m2)= ',
	1344	& qsurf(ig,igcm_h2o_ice,islope)
[38]	1345	write(,)' ==> Ice mesh North boundary (deg)= ',
	1346	& rlatv(j-1)*180./pi
[2909]	1347	enddo
[38]	1348	end if
	1349	enddo
	1350
	1351	c isotherm : Isothermal temperatures and no winds
	1352	c ------------------------------------------------
[164]	1353	else if (trim(modif) .eq. 'isotherm') then
[38]	1354
	1355	write(,)'Isothermal temperature of the atmosphere,
	1356	& surface and subsurface'
	1357	write(,) 'Value of this temperature ? :'
	1358	203 read(,,iostat=ierr) Tiso
	1359	if(ierr.ne.0) goto 203
	1360
	1361	do ig=1, ngridmx
[2909]	1362	do islope=1,nslope
	1363	tsurf(ig,islope) = Tiso
	1364	enddo
[38]	1365	end do
	1366	do l=2,nsoilmx
	1367	do ig=1, ngridmx
[2909]	1368	do islope=1,nslope
	1369	tsoil(ig,l,islope) = Tiso
	1370	enddo
[38]	1371	end do
	1372	end do
	1373	flagiso=.true.
	1374	call initial0(llm*ip1jmp1,ucov)
	1375	call initial0(llm*ip1jm,vcov)
	1376	call initial0(ngridmx*(llm+1),q2)
	1377
	1378	c co2ice=0 : remove CO2 polar ice caps'
	1379	c ------------------------------------------------
[164]	1380	else if (trim(modif) .eq. 'co2ice=0') then
[38]	1381	do ig=1,ngridmx
[2909]	1382	do islope=1,nslope
	1383	qsurf(ig,igcm_co2,islope)=0
	1384	emis(ig,islope)=emis(ngridmx/2,islope)
	1385	enddo
[38]	1386	end do
	1387
	1388	! therm_ini_s: (re)-set soil thermal inertia to reference surface values
	1389	! ----------------------------------------------------------------------
	1390
[164]	1391	else if (trim(modif).eq.'therm_ini_s') then
[38]	1392	! write(,)"surfithfi(1):",surfithfi(1)
	1393	do isoil=1,nsoilmx
	1394	inertiedat(1:ngridmx,isoil)=surfithfi(1:ngridmx)
	1395	enddo
[2942]	1396	do islope = 1,nslope
	1397	inertiesoil(:,:,islope) = inertiedat(:,:)
	1398	enddo
[38]	1399	write(,)'OK: Soil thermal inertia has been reset to referenc
	1400	&e surface values'
	1401	! write(,)"inertiedat(1,1):",inertiedat(1,1)
	1402	ithfi(:,:)=inertiedat(:,:)
	1403	! recast ithfi() onto ith()
	1404	call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
	1405	! Check:
	1406	! do i=1,iip1
	1407	! do j=1,jjp1
	1408	! do isoil=1,nsoilmx
	1409	! write(77,*) i,j,isoil," ",ith(i,j,isoil)
	1410	! enddo
	1411	! enddo
	1412	! enddo
	1413
	1414	! subsoilice_n: Put deep ice layer in northern hemisphere soil
	1415	! ------------------------------------------------------------
	1416
[164]	1417	else if (trim(modif).eq.'subsoilice_n') then
[38]	1418
	1419	write(,)'From which latitude (in deg.), up to the north pole,
	1420	&should we put subterranean ice?'
	1421	ierr=1
	1422	do while (ierr.ne.0)
	1423	read(,,iostat=ierr) val
	1424	if (ierr.eq.0) then ! got a value
	1425	! do a sanity check
	1426	if((val.lt.0.).or.(val.gt.90)) then
	1427	write(,)'Latitude should be between 0 and 90 deg. !!!'
	1428	ierr=1
	1429	else ! find corresponding jref (nearest latitude)
	1430	! note: rlatu(:) contains decreasing values of latitude
	1431	! starting from PI/2 to -PI/2
	1432	do j=1,jjp1
	1433	if ((rlatu(j)*180./pi.ge.val).and.
	1434	& (rlatu(j+1)*180./pi.le.val)) then
	1435	! find which grid point is nearest to val:
	1436	if (abs(rlatu(j)*180./pi-val).le.
	1437	& abs((rlatu(j+1)*180./pi-val))) then
	1438	jref=j
	1439	else
	1440	jref=j+1
	1441	endif
	1442
	1443	write(,)'Will use nearest grid latitude which is:',
	1444	& rlatu(jref)*180./pi
	1445	endif
	1446	enddo ! of do j=1,jjp1
	1447	endif ! of if((val.lt.0.).or.(val.gt.90))
	1448	endif !of if (ierr.eq.0)
	1449	enddo ! of do while
	1450
	1451	! Build layers() (as in soil_settings.F)
	1452	val2=sqrt(mlayer(0)*mlayer(1))
	1453	val3=mlayer(1)/mlayer(0)
	1454	do isoil=1,nsoilmx
	1455	layer(isoil)=val2(val3*(isoil-1))
	1456	enddo
	1457
	1458	write(,)'At which depth (in m.) does the ice layer begin?'
	1459	write(,)'(currently, the deepest soil layer extends down to:'
	1460	& ,layer(nsoilmx),')'
	1461	ierr=1
	1462	do while (ierr.ne.0)
	1463	read(,,iostat=ierr) val2
	1464	! write(,)'val2:',val2,'ierr=',ierr
	1465	if (ierr.eq.0) then ! got a value, but do a sanity check
	1466	if(val2.gt.layer(nsoilmx)) then
	1467	write(,)'Depth should be less than ',layer(nsoilmx)
	1468	ierr=1
	1469	endif
	1470	if(val2.lt.layer(1)) then
	1471	write(,)'Depth should be more than ',layer(1)
	1472	ierr=1
	1473	endif
	1474	endif
	1475	enddo ! of do while
	1476
	1477	! find the reference index iref the depth corresponds to
	1478	! if (val2.lt.layer(1)) then
	1479	! iref=1
	1480	! else
	1481	do isoil=1,nsoilmx-1
	1482	if((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1)))
	1483	& then
	1484	iref=isoil
	1485	exit
	1486	endif
	1487	enddo
	1488	! endif
	1489
	1490	! write(,)'iref:',iref,' jref:',jref
	1491	! write(,)'layer',layer
	1492	! write(,)'mlayer',mlayer
	1493
	1494	! thermal inertia of the ice:
	1495	ierr=1
	1496	do while (ierr.ne.0)
	1497	write(,)'What is the value of subterranean ice thermal inert
	1498	&ia? (e.g.: 2000)'
	1499	read(,,iostat=ierr)iceith
	1500	enddo ! of do while
	1501
	1502	! recast ithfi() onto ith()
	1503	call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
	1504
	1505	do j=1,jref
	1506	! write(,)'j:',j,'rlatu(j)180./pi:',rlatu(j)180./pi
	1507	do i=1,iip1 ! loop on longitudes
	1508	! Build "equivalent" thermal inertia for the mixed layer
	1509	ith(i,j,iref+1)=sqrt((layer(iref+1)-layer(iref))/
	1510	& (((val2-layer(iref))/(ith(i,j,iref)**2))+
	1511	& ((layer(iref+1)-val2)/(iceith)**2)))
	1512	! Set thermal inertia of lower layers
	1513	do isoil=iref+2,nsoilmx
	1514	ith(i,j,isoil)=iceith ! ice
	1515	enddo
	1516	enddo ! of do i=1,iip1
	1517	enddo ! of do j=1,jjp1
	1518
	1519
	1520	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
	1521
	1522	! do i=1,nsoilmx
	1523	! write(,)'i:',i,'ithfi(1,i):',ithfi(1,i)
	1524	! enddo
	1525
	1526
	1527	! subsoilice_s: Put deep ice layer in southern hemisphere soil
	1528	! ------------------------------------------------------------
	1529
[164]	1530	else if (trim(modif).eq.'subsoilice_s') then
[38]	1531
	1532	write(,)'From which latitude (in deg.), down to the south pol
	1533	&e, should we put subterranean ice?'
	1534	ierr=1
	1535	do while (ierr.ne.0)
	1536	read(,,iostat=ierr) val
	1537	if (ierr.eq.0) then ! got a value
	1538	! do a sanity check
	1539	if((val.gt.0.).or.(val.lt.-90)) then
	1540	write(,)'Latitude should be between 0 and -90 deg. !!!'
	1541	ierr=1
	1542	else ! find corresponding jref (nearest latitude)
	1543	! note: rlatu(:) contains decreasing values of latitude
	1544	! starting from PI/2 to -PI/2
	1545	do j=1,jjp1
	1546	if ((rlatu(j)*180./pi.ge.val).and.
	1547	& (rlatu(j+1)*180./pi.le.val)) then
	1548	! find which grid point is nearest to val:
	1549	if (abs(rlatu(j)*180./pi-val).le.
	1550	& abs((rlatu(j+1)*180./pi-val))) then
	1551	jref=j
	1552	else
	1553	jref=j+1
	1554	endif
	1555
	1556	write(,)'Will use nearest grid latitude which is:',
	1557	& rlatu(jref)*180./pi
	1558	endif
	1559	enddo ! of do j=1,jjp1
	1560	endif ! of if((val.lt.0.).or.(val.gt.90))
	1561	endif !of if (ierr.eq.0)
	1562	enddo ! of do while
	1563
	1564	! Build layers() (as in soil_settings.F)
	1565	val2=sqrt(mlayer(0)*mlayer(1))
	1566	val3=mlayer(1)/mlayer(0)
	1567	do isoil=1,nsoilmx
	1568	layer(isoil)=val2(val3*(isoil-1))
	1569	enddo
	1570
	1571	write(,)'At which depth (in m.) does the ice layer begin?'
	1572	write(,)'(currently, the deepest soil layer extends down to:'
	1573	& ,layer(nsoilmx),')'
	1574	ierr=1
	1575	do while (ierr.ne.0)
	1576	read(,,iostat=ierr) val2
	1577	! write(,)'val2:',val2,'ierr=',ierr
	1578	if (ierr.eq.0) then ! got a value, but do a sanity check
	1579	if(val2.gt.layer(nsoilmx)) then
	1580	write(,)'Depth should be less than ',layer(nsoilmx)
	1581	ierr=1
	1582	endif
	1583	if(val2.lt.layer(1)) then
	1584	write(,)'Depth should be more than ',layer(1)
	1585	ierr=1
	1586	endif
	1587	endif
	1588	enddo ! of do while
	1589
	1590	! find the reference index iref the depth corresponds to
	1591	do isoil=1,nsoilmx-1
	1592	if((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1)))
	1593	& then
	1594	iref=isoil
	1595	exit
	1596	endif
	1597	enddo
	1598
	1599	! write(,)'iref:',iref,' jref:',jref
	1600
	1601	! thermal inertia of the ice:
	1602	ierr=1
	1603	do while (ierr.ne.0)
	1604	write(,)'What is the value of subterranean ice thermal inert
	1605	&ia? (e.g.: 2000)'
	1606	read(,,iostat=ierr)iceith
	1607	enddo ! of do while
	1608
	1609	! recast ithfi() onto ith()
	1610	call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
	1611
	1612	do j=jref,jjp1
	1613	! write(,)'j:',j,'rlatu(j)180./pi:',rlatu(j)180./pi
	1614	do i=1,iip1 ! loop on longitudes
	1615	! Build "equivalent" thermal inertia for the mixed layer
	1616	ith(i,j,iref+1)=sqrt((layer(iref+1)-layer(iref))/
	1617	& (((val2-layer(iref))/(ith(i,j,iref)**2))+
	1618	& ((layer(iref+1)-val2)/(iceith)**2)))
	1619	! Set thermal inertia of lower layers
	1620	do isoil=iref+2,nsoilmx
	1621	ith(i,j,isoil)=iceith ! ice
	1622	enddo
	1623	enddo ! of do i=1,iip1
	1624	enddo ! of do j=jref,jjp1
	1625
	1626
	1627	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
	1628
	1629	c 'mons_ice' : use MONS data to build subsurface ice table
	1630	c --------------------------------------------------------
[164]	1631	else if (trim(modif).eq.'mons_ice') then
[38]	1632
	1633	! 1. Load MONS data
	1634	call load_MONS_data(MONS_Hdn,MONS_d21)
	1635
	1636	! 2. Get parameters from user
	1637	ierr=1
	1638	do while (ierr.ne.0)
	1639	write(,) "Coefficient to apply to MONS 'depth' in Northern",
	1640	& " Hemisphere?"
	1641	write(,) " (should be somewhere between 3.2e-4 and 1.3e-3)"
	1642	read(,,iostat=ierr) MONS_coeffN
	1643	enddo
	1644	ierr=1
	1645	do while (ierr.ne.0)
	1646	write(,) "Coefficient to apply to MONS 'depth' in Southern",
	1647	& " Hemisphere?"
	1648	write(,) " (should be somewhere between 3.2e-4 and 1.3e-3)"
	1649	read(,,iostat=ierr) MONS_coeffS
	1650	enddo
	1651	ierr=1
	1652	do while (ierr.ne.0)
	1653	write(,) "Value of subterranean ice thermal inertia ",
	1654	& " in Northern hemisphere?"
	1655	write(,) " (e.g.: 2000, or perhaps 2290)"
	1656	! read(,,iostat=ierr) iceith
	1657	read(,,iostat=ierr) iceithN
	1658	enddo
	1659	ierr=1
	1660	do while (ierr.ne.0)
	1661	write(,) "Value of subterranean ice thermal inertia ",
	1662	& " in Southern hemisphere?"
	1663	write(,) " (e.g.: 2000, or perhaps 2290)"
	1664	! read(,,iostat=ierr) iceith
	1665	read(,,iostat=ierr) iceithS
	1666	enddo
	1667
	1668	! 3. Build subterranean thermal inertia
	1669
	1670	! initialise subsurface inertia with reference surface values
	1671	do isoil=1,nsoilmx
	1672	ithfi(1:ngridmx,isoil)=surfithfi(1:ngridmx)
	1673	enddo
	1674	! recast ithfi() onto ith()
	1675	call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
	1676
	1677	do i=1,iip1 ! loop on longitudes
	1678	do j=1,jjp1 ! loop on latitudes
	1679	! set MONS_coeff
	1680	if (rlatu(j).ge.0) then ! northern hemisphere
	1681	! N.B: rlatu(:) contains decreasing values of latitude
	1682	! starting from PI/2 to -PI/2
	1683	MONS_coeff=MONS_coeffN
	1684	iceith=iceithN
	1685	else ! southern hemisphere
	1686	MONS_coeff=MONS_coeffS
	1687	iceith=iceithS
	1688	endif
	1689	! check if we should put subterranean ice
	1690	if (MONS_Hdn(i,j).ge.14.0) then ! no ice if Hdn<14%
	1691	! compute depth at which ice lies:
	1692	val=MONS_d21(i,j)*MONS_coeff
	1693	! compute val2= the diurnal skin depth of surface inertia
	1694	! assuming a volumetric heat cap. of C=1.e6 J.m-3.K-1
	1695	val2=ith(i,j,1)1.e-6sqrt(88775./3.14159)
	1696	if (val.lt.val2) then
	1697	! ice must be below the (surface inertia) diurnal skin depth
	1698	val=val2
	1699	endif
	1700	if (val.lt.layer(nsoilmx)) then ! subterranean ice
	1701	! find the reference index iref that depth corresponds to
	1702	iref=0
	1703	do isoil=1,nsoilmx-1
	1704	if ((val.ge.layer(isoil)).and.(val.lt.layer(isoil+1)))
	1705	& then
	1706	iref=isoil
	1707	exit
	1708	endif
	1709	enddo
	1710	! Build "equivalent" thermal inertia for the mixed layer
	1711	ith(i,j,iref+1)=sqrt((layer(iref+1)-layer(iref))/
	1712	& (((val-layer(iref))/(ith(i,j,iref+1)**2))+
	1713	& ((layer(iref+1)-val)/(iceith)**2)))
	1714	! Set thermal inertia of lower layers
	1715	do isoil=iref+2,nsoilmx
	1716	ith(i,j,isoil)=iceith
	1717	enddo
	1718	endif ! of if (val.lt.layer(nsoilmx))
	1719	endif ! of if (MONS_Hdn(i,j).lt.14.0)
	1720	enddo ! do j=1,jjp1
	1721	enddo ! do i=1,iip1
	1722
	1723	! Check:
	1724	! do i=1,iip1
	1725	! do j=1,jjp1
	1726	! do isoil=1,nsoilmx
	1727	! write(77,*) i,j,isoil," ",ith(i,j,isoil)
	1728	! enddo
	1729	! enddo
	1730	! enddo
	1731
	1732	! recast ith() into ithfi()
	1733	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
[2909]	1734
	1735	else if (trim(modif) .eq. 'nslope') then
	1736	write(,) 'set a new number of subgrid scale slope'
	1737	write(,) 'Current value=', nslope
	1738	write(,) 'Enter value for nslope (ex: 1,5,7)?'
	1739	ierr=1
	1740	do while (ierr.ne.0)
	1741	read(,,iostat=ierr) nslope_new
	1742	enddo
	1743
	1744	write(,) 'This can take some time...'
	1745	write(,) 'You can go grab a coffee and relax a bit'
	1746
	1747	if(nslope.eq.nslope_new) then
	1748	write(,) 'The number of subslope you entered is the same'
	1749	write(,) 'as the number written in startfi.nc. '
	1750	write(,) 'Nothing will be done'
	1751	else
	1752
	1753	nslope_old=nslope
	1754	nslope=nslope_new
	1755
	1756	call end_comslope_h
	1757	call ini_comslope_h(ngridmx,nslope_new)
	1758
	1759	allocate(default_def_slope(nslope_new+1))
	1760	!Sub-grid scale subslopes
	1761	if (nslope_new.eq.7) then
	1762	default_def_slope(1) = -43.
	1763	default_def_slope(2) = -19.
	1764	default_def_slope(3) = -9.
	1765	default_def_slope(4) = -3.
	1766	default_def_slope(5) = 3.
	1767	default_def_slope(6) = 9.
	1768	default_def_slope(7) = 19.
	1769	default_def_slope(8) = 43.
	1770	elseif (nslope_new.eq.5) then
	1771	default_def_slope(1) = -43.
	1772	default_def_slope(2) = -9.
	1773	default_def_slope(3) = -3.
	1774	default_def_slope(4) = 3.
	1775	default_def_slope(5) = 9.
	1776	default_def_slope(6) = 43.
	1777	elseif (nslope_new.eq.1) then
	1778	default_def_slope(1) = -50.
	1779	default_def_slope(2) = 50.
	1780	endif
	1781
	1782	do islope=1,nslope_new+1
	1783	def_slope(islope) = default_def_slope(islope)
	1784	enddo
	1785
	1786	do islope=1,nslope_new
	1787	def_slope_mean(islope) =(def_slope(islope)+def_slope(islope+1))/2.
	1788	enddo
	1789
	1790	iflat = 1
	1791	DO islope=2,nslope_new
	1792	IF(abs(def_slope_mean(islope)).lt.
	1793	& abs(def_slope_mean(iflat)))THEN
	1794	iflat = islope
	1795	ENDIF
	1796	ENDDO
	1797
	1798	call subslope_mola(ngridmx,nslope_new,def_slope,subslope_dist)
	1799
	1800	! Surfdat related stuff
[38]	1801
[2909]	1802	allocate(tsurf_old_slope(ngridmx,nslope_old))
	1803	allocate(qsurf_old_slope(ngridmx,nqtot,nslope_old))
	1804	allocate(emis_old_slope(ngridmx,nslope_old))
	1805	allocate(watercap_old_slope(ngridmx,nslope_old))
[2999]	1806	allocate(perenial_co2_old_slope(ngridmx,nslope_old))
[2909]	1807
[2952]	1808
[2909]	1809	tsurf_old_slope(:,:)=tsurf(:,:)
	1810	qsurf_old_slope(:,:,:)=qsurf(:,:,:)
	1811	emis_old_slope(:,:)=emis(:,:)
	1812	watercap_old_slope(:,:)=watercap(:,:)
[2999]	1813	perenial_co2_old_slope(:,:) = perenial_co2ice(:,:)
[2909]	1814	call end_surfdat_h_slope_var
	1815	call ini_surfdat_h_slope_var(ngridmx,nqtot,nslope_new)
	1816
	1817	! Comsoil related stuff (tsoil)
	1818
	1819	allocate(tsoil_old_slope(ngridmx,nsoilmx,nslope_old))
[2952]	1820	allocate(inertiesoil_old_slope(ngridmx,nsoilmx,nslope_old))
	1821	allocate(flux_geo_old_slope(ngridmx,nslope_old))
[2909]	1822
[2952]	1823	inertiesoil_old_slope(:,:,:)=inertiesoil(:,:,:)
[2909]	1824	tsoil_old_slope(:,:,:)=tsoil(:,:,:)
[2952]	1825	flux_geo_old_slope(:,:)=flux_geo(:,:)
[2909]	1826
	1827	call end_comsoil_h_slope_var
	1828	call ini_comsoil_h_slope_var(ngridmx,nslope_new)
	1829
	1830	! Dimradmars related stuff (albedo)
	1831
	1832	allocate(albedo_old_slope(ngridmx,2,nslope_old))
	1833	albedo_old_slope(:,:,:)=albedo(:,:,:)
	1834
	1835	call end_dimradmars_mod_slope_var
	1836	call ini_dimradmars_mod_slope_var(ngridmx,nslope_new)
	1837
	1838	if(nslope_old.eq.1 .and. nslope_new.gt.1) then
	1839	do islope=1,nslope_new
	1840	tsurf(:,islope)=tsurf_old_slope(:,1)
	1841	qsurf(:,:,islope)=qsurf_old_slope(:,:,1)
	1842	emis(:,islope)=emis_old_slope(:,1)
	1843	watercap(:,islope)=watercap_old_slope(:,1)
[2999]	1844	perenial_co2ice(:,islope)= perenial_co2_old_slope(:,1)
[2909]	1845	tsoil(:,:,islope)=tsoil_old_slope(:,:,1)
	1846	albedo(:,:,islope)=albedo_old_slope(:,:,1)
[2952]	1847	inertiesoil(:,:,islope)=inertiesoil_old_slope(:,:,1)
	1848	flux_geo(:,islope)=flux_geo_old_slope(:,1)
[2909]	1849	enddo
	1850	elseif(nslope_new.eq.1) then
	1851	tsurf(:,1)=tsurf_old_slope(:,iflat)
	1852	qsurf(:,:,1)=qsurf_old_slope(:,:,iflat)
	1853	emis(:,1)=emis_old_slope(:,iflat)
	1854	watercap(:,1)=watercap_old_slope(:,iflat)
[2999]	1855	perenial_co2ice(:,islope)= perenial_co2_old_slope(:,iflat)
[2909]	1856	tsoil(:,:,1)=tsoil_old_slope(:,:,iflat)
	1857	albedo(:,:,1)=albedo_old_slope(:,:,iflat)
[2952]	1858	inertiesoil(:,:,1)=inertiesoil_old_slope(:,:,iflat)
	1859	flux_geo(:,1)=flux_geo_old_slope(:,iflat)
[2909]	1860	elseif(nslope_old.eq.5 .and. nslope_new.eq.7) then
	1861	do islope=1,nslope_new
	1862	tsurf(:,islope)=tsurf_old_slope(:,iflat)
	1863	qsurf(:,:,islope)=qsurf_old_slope(:,:,iflat)
	1864	emis(:,islope)=emis_old_slope(:,iflat)
	1865	watercap(:,islope)=watercap_old_slope(:,iflat)
[2999]	1866	perenial_co2ice(:,islope)= perenial_co2_old_slope(:,iflat)
[2909]	1867	tsoil(:,:,islope)=tsoil_old_slope(:,:,iflat)
	1868	albedo(:,:,islope)=albedo_old_slope(:,:,iflat)
[2952]	1869	inertiesoil(:,:,islope)=inertiesoil_old_slope(:,:,iflat)
	1870	flux_geo(:,islope)=flux_geo_old_slope(:,iflat)
[2909]	1871	enddo
	1872	elseif(nslope_old.eq.7 .and. nslope_new.eq.5) then
	1873	do islope=1,nslope_new
	1874	tsurf(:,islope)=tsurf_old_slope(:,iflat)
	1875	qsurf(:,:,islope)=qsurf_old_slope(:,:,iflat)
	1876	emis(:,islope)=emis_old_slope(:,iflat)
	1877	watercap(:,islope)=watercap_old_slope(:,iflat)
[2999]	1878	perenial_co2ice(:,islope)= perenial_co2_old_slope(:,iflat)
[2909]	1879	tsoil(:,:,islope)=tsoil_old_slope(:,:,iflat)
	1880	albedo(:,:,islope)=albedo_old_slope(:,:,iflat)
[2952]	1881	inertiesoil(:,:,islope)=inertiesoil_old_slope(:,:,iflat)
	1882	flux_geo(:,islope)=flux_geo_old_slope(:,iflat)
[2909]	1883	enddo
	1884	else
	1885	write(,)' Problem choice of nslope'
	1886	write(,)' Value not taken into account'
	1887	CALL ABORT
	1888	endif
	1889
	1890	endif !nslope=nslope_new
	1891
[38]	1892	else
	1893	write(,) ' Unknown (misspelled?) option!!!'
[164]	1894	end if ! of if (trim(modif) .eq. '...') elseif ...
[38]	1895
	1896	enddo ! of do ! infinite loop on liste of changes
	1897
	1898	999 continue
	1899
	1900
	1901	c=======================================================================
	1902	c Correct pressure on the new grid (menu 0)
	1903	c=======================================================================
	1904
	1905	if (choix_1.eq.0) then
	1906	r = 1000.*8.31/mugaz
	1907
	1908	do j=1,jjp1
	1909	do i=1,iip1
	1910	ps(i,j) = ps(i,j) *
	1911	. exp((phisold_newgrid(i,j)-phis(i,j)) /
	1912	. (t(i,j,1) * r))
	1913	end do
	1914	end do
	1915
	1916	c periodicity of surface ps in longitude
	1917	do j=1,jjp1
	1918	ps(1,j) = ps(iip1,j)
	1919	end do
	1920	end if
	1921
	1922	c=======================================================================
	1923	c=======================================================================
	1924
	1925	c=======================================================================
	1926	c Initialisation de la physique / ecriture de newstartfi :
	1927	c=======================================================================
	1928
	1929
	1930	CALL inifilr
	1931	CALL pression(ip1jmp1, ap, bp, ps, p3d)
	1932
	1933	c-----------------------------------------------------------------------
	1934	c Initialisation pks:
	1935	c-----------------------------------------------------------------------
	1936
[2167]	1937	CALL exner_hyb(ip1jmp1, ps, p3d, pks, pk, pkf)
[38]	1938	! Calcul de la temperature potentielle teta
	1939
	1940	if (flagiso) then
	1941	DO l=1,llm
	1942	DO j=1,jjp1
	1943	DO i=1,iim
	1944	teta(i,j,l) = Tiso * cpp/pk(i,j,l)
	1945	ENDDO
	1946	teta (iip1,j,l)= teta (1,j,l)
	1947	ENDDO
	1948	ENDDO
	1949	else if (choix_1.eq.0) then
	1950	DO l=1,llm
	1951	DO j=1,jjp1
	1952	DO i=1,iim
	1953	teta(i,j,l) = t(i,j,l) * cpp/pk(i,j,l)
	1954	ENDDO
	1955	teta (iip1,j,l)= teta (1,j,l)
	1956	ENDDO
	1957	ENDDO
	1958	endif
	1959
	1960	C Calcul intermediaire
	1961	c
	1962	if (choix_1.eq.0) then
	1963	CALL massdair( p3d, masse )
	1964	c
	1965	print *,' ALPHAX ',alphax
	1966
	1967	DO l = 1, llm
	1968	DO i = 1, iim
	1969	xppn(i) = aire( i, 1 ) * masse( i , 1 , l )
	1970	xpps(i) = aire( i,jjp1 ) * masse( i , jjp1 , l )
	1971	ENDDO
	1972	xpn = SUM(xppn)/apoln
	1973	xps = SUM(xpps)/apols
	1974	DO i = 1, iip1
	1975	masse( i , 1 , l ) = xpn
	1976	masse( i , jjp1 , l ) = xps
	1977	ENDDO
	1978	ENDDO
	1979	endif
	1980	phis(iip1,:) = phis(1,:)
	1981
[1415]	1982	c CALL inidissip ( lstardis, nitergdiv, nitergrot, niterh,
	1983	c * tetagdiv, tetagrot , tetatemp )
[38]	1984	itau=0
	1985	if (choix_1.eq.0) then
	1986	day_ini=int(date)
[999]	1987	hour_ini=date-int(date)
[38]	1988	endif
	1989	c
	1990	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
	1991
	1992	CALL caldyn0( itau,ucov,vcov,teta,ps,masse,pk,phis ,
	1993	* phi,w, pbaru,pbarv,day_ini+time )
	1994	c CALL caldyn
	1995	c $ ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
	1996	c $ phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, day_ini )
	1997
[1415]	1998	CALL dynredem0("restart.nc",day_ini,phis)
[999]	1999	CALL dynredem1("restart.nc",hour_ini,vcov,ucov,teta,q,
[1415]	2000	. masse,ps)
[38]	2001	C
	2002	C Ecriture etat initial physique
	2003	C
[2001]	2004	call physdem0("restartfi.nc",longitude,latitude,
	2005	& nsoilmx,ngridmx,llm,
	2006	& nqtot,dtphys,real(day_ini),0.0,cell_area,
[2909]	2007	& albfi,ithfi,def_slope,
	2008	& subslope_dist)
[1047]	2009	call physdem1("restartfi.nc",nsoilmx,ngridmx,llm,nqtot,
[1944]	2010	& dtphys,hour_ini,
[2943]	2011	& tsurf,tsoil,inertiesoil,albedo,emis,q2,qsurf,
[2999]	2012	& tauscaling,totcloudfrac,wstar,watercap,
	2013	& perenial_co2ice)
[38]	2014
	2015	c=======================================================================
	2016	c Formats
	2017	c=======================================================================
	2018
	2019	1 FORMAT(//10x,'la valeur de im =',i4,2x,'lue sur le fichier de dema
	2020	*rrage est differente de la valeur parametree iim =',i4//)
	2021	2 FORMAT(//10x,'la valeur de jm =',i4,2x,'lue sur le fichier de dema
	2022	*rrage est differente de la valeur parametree jjm =',i4//)
	2023	3 FORMAT(//10x,'la valeur de lllm =',i4,2x,'lue sur le fichier demar
	2024	*rage est differente de la valeur parametree llm =',i4//)
	2025
[171]	2026	write(,) "newstart: All is well that ends well."
	2027
[38]	2028	end
	2029
	2030	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2031	subroutine load_MONS_data(MONS_Hdn,MONS_d21)
[1918]	2032
	2033	use datafile_mod, only:datadir
	2034
[38]	2035	implicit none
	2036	! routine to load Benedicte Diez MONS dataset, fill in date in southern
	2037	! polar region, and interpolate the result onto the GCM grid
[1918]	2038	include"dimensions.h"
	2039	include"paramet.h"
	2040	include"comgeom.h"
[38]	2041	! arguments:
	2042	real,intent(out) :: MONS_Hdn(iip1,jjp1) ! Hdn: %WEH=Mass fraction of H2O
	2043	real,intent(out) :: MONS_d21(iip1,jjp1) ! ice table "depth" (in kg/m2)
	2044	! N.B MONS datasets should be of dimension (iip1,jjp1)
	2045	! local variables:
	2046	character(len=88) :: filename="results_MONS_lat_lon_H_depth.txt"
	2047	character(len=88) :: txt ! to store some text
	2048	integer :: ierr,i,j
	2049	integer,parameter :: nblon=180 ! number of longitudes of MONS datasets
	2050	integer,parameter :: nblat=90 ! number of latitudes of MONS datasets
	2051	real :: pi
	2052	real :: longitudes(nblon) ! MONS dataset longitudes
	2053	real :: latitudes(nblat) ! MONS dataset latitudes
	2054	! MONS dataset: mass fraction of H2O where H is assumed to be in H2O
	2055	real :: Hdn(nblon,nblat)
	2056	real :: d21(nblon,nblat)! MONS dataset "depth" (g/cm2)
	2057
	2058	! Extended MONS dataset (for interp_horiz)
	2059	real :: Hdnx(nblon+1,nblat)
	2060	real :: d21x(nblon+1,nblat)
	2061	real :: lon_bound(nblon+1) ! longitude boundaries
	2062	real :: lat_bound(nblat-1) ! latitude boundaries
	2063
	2064	! 1. Initializations:
	2065
	2066	write(,) "Loading MONS data"
	2067
	2068	! Open MONS datafile:
[1918]	2069	open(42,file=trim(datadir)//"/"//trim(filename),
[38]	2070	& status="old",iostat=ierr)
	2071	if (ierr/=0) then
	2072	write(,) "Error in load_MONS_data:"
	2073	write(,) "Failed opening file ",
[1918]	2074	& trim(datadir)//"/"//trim(filename)
	2075	write(,)'1) You can change this directory address in ',
	2076	& 'callfis.def with'
	2077	write(,)' datadir=/path/to/datafiles'
[38]	2078	write(,)'2) If necessary ',trim(filename),
	2079	& ' (and other datafiles)'
	2080	write(,)' can be obtained online at:'
[1383]	2081	write(,)'http://www.lmd.jussieu.fr/~lmdz/planets/mars/datadir'
[38]	2082	CALL ABORT
	2083	else ! skip first line of file (dummy read)
	2084	read(42,*) txt
	2085	endif
	2086
	2087	pi=2.*asin(1.)
	2088
	2089	!2. Load MONS data (on MONS grid)
	2090	do j=1,nblat
	2091	do i=1,nblon
	2092	! swap latitude index so latitudes go from north pole to south pole:
	2093	read(42,*) latitudes(nblat-j+1),longitudes(i),
	2094	& Hdn(i,nblat-j+1),d21(i,nblat-j+1)
	2095	! multiply d21 by 10 to convert from g/cm2 to kg/m2
	2096	d21(i,nblat-j+1)=d21(i,nblat-j+1)*10.0
	2097	enddo
	2098	enddo
	2099	close(42)
	2100
	2101	! there is unfortunately no d21 data for latitudes -77 to -90
	2102	! so we build some by linear interpolation between values at -75
	2103	! and assuming d21=0 at the pole
	2104	do j=84,90 ! latitudes(84)=-77 ; latitudes(83)=-75
	2105	do i=1,nblon
	2106	d21(i,j)=d21(i,83)*((latitudes(j)+90)/15.0)
	2107	enddo
	2108	enddo
	2109
	2110	! 3. Build extended MONS dataset & boundaries (for interp_horiz)
	2111	! longitude boundaries (in radians):
	2112	do i=1,nblon
	2113	! NB: MONS data is every 2 degrees in longitude
	2114	lon_bound(i)=(longitudes(i)+1.0)*pi/180.0
	2115	enddo
	2116	! extra 'modulo' value
	2117	lon_bound(nblon+1)=lon_bound(1)+2.0*pi
	2118
	2119	! latitude boundaries (in radians):
	2120	do j=1,nblat-1
	2121	! NB: Mons data is every 2 degrees in latitude
	2122	lat_bound(j)=(latitudes(j)-1.0)*pi/180.0
	2123	enddo
	2124
	2125	! MONS datasets:
	2126	do j=1,nblat
	2127	Hdnx(1:nblon,j)=Hdn(1:nblon,j)
	2128	Hdnx(nblon+1,j)=Hdnx(1,j)
	2129	d21x(1:nblon,j)=d21(1:nblon,j)
	2130	d21x(nblon+1,j)=d21x(1,j)
	2131	enddo
	2132
	2133	! Interpolate onto GCM grid
	2134	call interp_horiz(Hdnx,MONS_Hdn,nblon,nblat-1,iim,jjm,1,
	2135	& lon_bound,lat_bound,rlonu,rlatv)
	2136	call interp_horiz(d21x,MONS_d21,nblon,nblat-1,iim,jjm,1,
	2137	& lon_bound,lat_bound,rlonu,rlatv)
	2138
	2139	end subroutine

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

Download in other formats:

Original Format