Context Navigation

newstart.F @ 3571

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

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

Download in other formats:

Original Format