Context Navigation

newstart.F @ 3302

Last change on this file since 3302 was 3203, checked in by jbclement, 11 months ago

Mars PCM:

Addition of the possibility to use subslopes parametization in 1D. To do so, put 'nslope=x' in "run.def" where x (1, 3, 5 or 7) is the number of slopes you want to, or modify it in an appropriate "startfi.nc". Then, default subslopes definition and distribution are used by the model. The already existing case of using one slope whose inclination and orientation are defined in "run.def" is still possible.
Some cleanings throughout the code, in particular related to unused variables.

JBC

File size: 76.8 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
[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)
[3139]	200	REAL,ALLOCATABLE :: perennial_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,
[3139]	447	& tauscaling,totcloudfrac,surfith,nid,watercap,perennial_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,
[3118]	467	& nqsoil,day_ini,time,tsurf,tsoil,albedo,emis,
	468	& q2,qsurf,qsoil,tauscaling,totcloudfrac,
[3139]	469	& wstar,watercap,perennial_co2ice,
[2999]	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
[3183]	1736	write(,) 'set a new number of subgrid scale slopes'
[2909]	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))
[3183]	1760	! Sub-grid scale slopes parameters (minimum/maximun angles)
	1761	select case(nslope)
	1762	case(1)
[3203]	1763	default_def_slope(1) = -90.
	1764	default_def_slope(2) = 90.
	1765	case(3)
[3183]	1766	default_def_slope(1) = -50.
[3203]	1767	default_def_slope(2) = -3.
	1768	default_def_slope(3) = 3.
	1769	default_def_slope(4) = 50.
[3183]	1770	case(5)
	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	case(7)
	1778	default_def_slope(1) = -43.
	1779	default_def_slope(2) = -19.
	1780	default_def_slope(3) = -9.
	1781	default_def_slope(4) = -3.
	1782	default_def_slope(5) = 3.
	1783	default_def_slope(6) = 9.
	1784	default_def_slope(7) = 19.
	1785	default_def_slope(8) = 43.
	1786	case default
	1787	write(,) 'Number of slopes not possible: nslope should
[3203]	1788	& be 1, 3, 5 or 7!'
[3183]	1789	call abort
	1790	end select
[2909]	1791
	1792	do islope=1,nslope_new+1
	1793	def_slope(islope) = default_def_slope(islope)
	1794	enddo
	1795
	1796	do islope=1,nslope_new
[3183]	1797	def_slope_mean(islope)=(def_slope(islope)+def_slope(islope+1))/2.
[2909]	1798	enddo
	1799
	1800	iflat = 1
	1801	DO islope=2,nslope_new
	1802	IF(abs(def_slope_mean(islope)).lt.
	1803	& abs(def_slope_mean(iflat)))THEN
	1804	iflat = islope
	1805	ENDIF
	1806	ENDDO
	1807
[3203]	1808	if (ngridmx /= 1) then
	1809	call subslope_mola(ngridmx,nslope_new,def_slope,subslope_dist)
	1810	endif
[2909]	1811
[3183]	1812	! Surfdat related stuff
[2909]	1813	allocate(tsurf_old_slope(ngridmx,nslope_old))
	1814	allocate(qsurf_old_slope(ngridmx,nqtot,nslope_old))
	1815	allocate(emis_old_slope(ngridmx,nslope_old))
	1816	allocate(watercap_old_slope(ngridmx,nslope_old))
[3139]	1817	allocate(perennial_co2_old_slope(ngridmx,nslope_old))
[2909]	1818
	1819	tsurf_old_slope(:,:)=tsurf(:,:)
	1820	qsurf_old_slope(:,:,:)=qsurf(:,:,:)
	1821	emis_old_slope(:,:)=emis(:,:)
	1822	watercap_old_slope(:,:)=watercap(:,:)
[3139]	1823	perennial_co2_old_slope(:,:) = perennial_co2ice(:,:)
[2909]	1824	call end_surfdat_h_slope_var
	1825	call ini_surfdat_h_slope_var(ngridmx,nqtot,nslope_new)
	1826
	1827	! Comsoil related stuff (tsoil)
	1828	allocate(tsoil_old_slope(ngridmx,nsoilmx,nslope_old))
[2952]	1829	allocate(inertiesoil_old_slope(ngridmx,nsoilmx,nslope_old))
	1830	allocate(flux_geo_old_slope(ngridmx,nslope_old))
[2909]	1831
[2952]	1832	inertiesoil_old_slope(:,:,:)=inertiesoil(:,:,:)
[2909]	1833	tsoil_old_slope(:,:,:)=tsoil(:,:,:)
[2952]	1834	flux_geo_old_slope(:,:)=flux_geo(:,:)
[2909]	1835
	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
	1843	call end_dimradmars_mod_slope_var
	1844	call ini_dimradmars_mod_slope_var(ngridmx,nslope_new)
	1845
	1846	if(nslope_old.eq.1 .and. nslope_new.gt.1) then
	1847	do islope=1,nslope_new
	1848	tsurf(:,islope)=tsurf_old_slope(:,1)
	1849	qsurf(:,:,islope)=qsurf_old_slope(:,:,1)
	1850	emis(:,islope)=emis_old_slope(:,1)
	1851	watercap(:,islope)=watercap_old_slope(:,1)
[3139]	1852	perennial_co2ice(:,islope)= perennial_co2_old_slope(:,1)
[2909]	1853	tsoil(:,:,islope)=tsoil_old_slope(:,:,1)
	1854	albedo(:,:,islope)=albedo_old_slope(:,:,1)
[2952]	1855	inertiesoil(:,:,islope)=inertiesoil_old_slope(:,:,1)
	1856	flux_geo(:,islope)=flux_geo_old_slope(:,1)
[2909]	1857	enddo
	1858	elseif(nslope_new.eq.1) then
	1859	tsurf(:,1)=tsurf_old_slope(:,iflat)
	1860	qsurf(:,:,1)=qsurf_old_slope(:,:,iflat)
	1861	emis(:,1)=emis_old_slope(:,iflat)
	1862	watercap(:,1)=watercap_old_slope(:,iflat)
[3139]	1863	perennial_co2ice(:,islope)=
	1864	& perennial_co2_old_slope(:,iflat)
[2909]	1865	tsoil(:,:,1)=tsoil_old_slope(:,:,iflat)
	1866	albedo(:,:,1)=albedo_old_slope(:,:,iflat)
[2952]	1867	inertiesoil(:,:,1)=inertiesoil_old_slope(:,:,iflat)
	1868	flux_geo(:,1)=flux_geo_old_slope(:,iflat)
[2909]	1869	elseif(nslope_old.eq.5 .and. nslope_new.eq.7) then
	1870	do islope=1,nslope_new
	1871	tsurf(:,islope)=tsurf_old_slope(:,iflat)
	1872	qsurf(:,:,islope)=qsurf_old_slope(:,:,iflat)
	1873	emis(:,islope)=emis_old_slope(:,iflat)
	1874	watercap(:,islope)=watercap_old_slope(:,iflat)
[3139]	1875	perennial_co2ice(:,islope)=
	1876	& perennial_co2_old_slope(:,iflat)
[2909]	1877	tsoil(:,:,islope)=tsoil_old_slope(:,:,iflat)
	1878	albedo(:,:,islope)=albedo_old_slope(:,:,iflat)
[2952]	1879	inertiesoil(:,:,islope)=inertiesoil_old_slope(:,:,iflat)
	1880	flux_geo(:,islope)=flux_geo_old_slope(:,iflat)
[2909]	1881	enddo
	1882	elseif(nslope_old.eq.7 .and. nslope_new.eq.5) then
	1883	do islope=1,nslope_new
	1884	tsurf(:,islope)=tsurf_old_slope(:,iflat)
	1885	qsurf(:,:,islope)=qsurf_old_slope(:,:,iflat)
	1886	emis(:,islope)=emis_old_slope(:,iflat)
	1887	watercap(:,islope)=watercap_old_slope(:,iflat)
[3139]	1888	perennial_co2ice(:,islope)=
	1889	& perennial_co2_old_slope(:,iflat)
[2909]	1890	tsoil(:,:,islope)=tsoil_old_slope(:,:,iflat)
	1891	albedo(:,:,islope)=albedo_old_slope(:,:,iflat)
[2952]	1892	inertiesoil(:,:,islope)=inertiesoil_old_slope(:,:,iflat)
	1893	flux_geo(:,islope)=flux_geo_old_slope(:,iflat)
[2909]	1894	enddo
	1895	else
	1896	write(,)' Problem choice of nslope'
	1897	write(,)' Value not taken into account'
	1898	CALL ABORT
	1899	endif
	1900
	1901	endif !nslope=nslope_new
	1902
[38]	1903	else
	1904	write(,) ' Unknown (misspelled?) option!!!'
[164]	1905	end if ! of if (trim(modif) .eq. '...') elseif ...
[38]	1906
	1907	enddo ! of do ! infinite loop on liste of changes
	1908
	1909	999 continue
	1910
	1911
	1912	c=======================================================================
	1913	c Correct pressure on the new grid (menu 0)
	1914	c=======================================================================
	1915
	1916	if (choix_1.eq.0) then
	1917	r = 1000.*8.31/mugaz
	1918
	1919	do j=1,jjp1
	1920	do i=1,iip1
	1921	ps(i,j) = ps(i,j) *
	1922	. exp((phisold_newgrid(i,j)-phis(i,j)) /
	1923	. (t(i,j,1) * r))
	1924	end do
	1925	end do
	1926
	1927	c periodicity of surface ps in longitude
	1928	do j=1,jjp1
	1929	ps(1,j) = ps(iip1,j)
	1930	end do
	1931	end if
	1932
	1933	c=======================================================================
	1934	c=======================================================================
	1935
	1936	c=======================================================================
	1937	c Initialisation de la physique / ecriture de newstartfi :
	1938	c=======================================================================
	1939
	1940
	1941	CALL inifilr
	1942	CALL pression(ip1jmp1, ap, bp, ps, p3d)
	1943
	1944	c-----------------------------------------------------------------------
	1945	c Initialisation pks:
	1946	c-----------------------------------------------------------------------
	1947
[2167]	1948	CALL exner_hyb(ip1jmp1, ps, p3d, pks, pk, pkf)
[38]	1949	! Calcul de la temperature potentielle teta
	1950
	1951	if (flagiso) then
	1952	DO l=1,llm
	1953	DO j=1,jjp1
	1954	DO i=1,iim
	1955	teta(i,j,l) = Tiso * cpp/pk(i,j,l)
	1956	ENDDO
	1957	teta (iip1,j,l)= teta (1,j,l)
	1958	ENDDO
	1959	ENDDO
	1960	else if (choix_1.eq.0) then
	1961	DO l=1,llm
	1962	DO j=1,jjp1
	1963	DO i=1,iim
	1964	teta(i,j,l) = t(i,j,l) * cpp/pk(i,j,l)
	1965	ENDDO
	1966	teta (iip1,j,l)= teta (1,j,l)
	1967	ENDDO
	1968	ENDDO
	1969	endif
	1970
	1971	C Calcul intermediaire
	1972	c
	1973	if (choix_1.eq.0) then
	1974	CALL massdair( p3d, masse )
	1975	c
	1976	print *,' ALPHAX ',alphax
	1977
	1978	DO l = 1, llm
	1979	DO i = 1, iim
	1980	xppn(i) = aire( i, 1 ) * masse( i , 1 , l )
	1981	xpps(i) = aire( i,jjp1 ) * masse( i , jjp1 , l )
	1982	ENDDO
	1983	xpn = SUM(xppn)/apoln
	1984	xps = SUM(xpps)/apols
	1985	DO i = 1, iip1
	1986	masse( i , 1 , l ) = xpn
	1987	masse( i , jjp1 , l ) = xps
	1988	ENDDO
	1989	ENDDO
	1990	endif
	1991	phis(iip1,:) = phis(1,:)
	1992
[1415]	1993	c CALL inidissip ( lstardis, nitergdiv, nitergrot, niterh,
	1994	c * tetagdiv, tetagrot , tetatemp )
[38]	1995	itau=0
	1996	if (choix_1.eq.0) then
	1997	day_ini=int(date)
[999]	1998	hour_ini=date-int(date)
[38]	1999	endif
	2000	c
	2001	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
	2002
	2003	CALL caldyn0( itau,ucov,vcov,teta,ps,masse,pk,phis ,
	2004	* phi,w, pbaru,pbarv,day_ini+time )
	2005	c CALL caldyn
	2006	c $ ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
	2007	c $ phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, day_ini )
	2008
[1415]	2009	CALL dynredem0("restart.nc",day_ini,phis)
[999]	2010	CALL dynredem1("restart.nc",hour_ini,vcov,ucov,teta,q,
[1415]	2011	. masse,ps)
[38]	2012	C
	2013	C Ecriture etat initial physique
	2014	C
[2001]	2015	call physdem0("restartfi.nc",longitude,latitude,
	2016	& nsoilmx,ngridmx,llm,
	2017	& nqtot,dtphys,real(day_ini),0.0,cell_area,
[2909]	2018	& albfi,ithfi,def_slope,
	2019	& subslope_dist)
[3118]	2020	call physdem1("restartfi.nc",nsoilmx,ngridmx,llm,nqtot,nqsoil,
[1944]	2021	& dtphys,hour_ini,
[3118]	2022	& tsurf,tsoil,inertiesoil,albedo,
	2023	& emis,q2,qsurf,qsoil,
[2999]	2024	& tauscaling,totcloudfrac,wstar,watercap,
[3139]	2025	& perennial_co2ice)
[38]	2026
	2027	c=======================================================================
	2028	c Formats
	2029	c=======================================================================
	2030
	2031	1 FORMAT(//10x,'la valeur de im =',i4,2x,'lue sur le fichier de dema
	2032	*rrage est differente de la valeur parametree iim =',i4//)
	2033	2 FORMAT(//10x,'la valeur de jm =',i4,2x,'lue sur le fichier de dema
	2034	*rrage est differente de la valeur parametree jjm =',i4//)
	2035	3 FORMAT(//10x,'la valeur de lllm =',i4,2x,'lue sur le fichier demar
	2036	*rage est differente de la valeur parametree llm =',i4//)
	2037
[171]	2038	write(,) "newstart: All is well that ends well."
	2039
[38]	2040	end
	2041
	2042	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2043	subroutine load_MONS_data(MONS_Hdn,MONS_d21)
[1918]	2044
	2045	use datafile_mod, only:datadir
	2046
[38]	2047	implicit none
	2048	! routine to load Benedicte Diez MONS dataset, fill in date in southern
	2049	! polar region, and interpolate the result onto the GCM grid
[1918]	2050	include"dimensions.h"
	2051	include"paramet.h"
	2052	include"comgeom.h"
[38]	2053	! arguments:
	2054	real,intent(out) :: MONS_Hdn(iip1,jjp1) ! Hdn: %WEH=Mass fraction of H2O
	2055	real,intent(out) :: MONS_d21(iip1,jjp1) ! ice table "depth" (in kg/m2)
	2056	! N.B MONS datasets should be of dimension (iip1,jjp1)
	2057	! local variables:
	2058	character(len=88) :: filename="results_MONS_lat_lon_H_depth.txt"
	2059	character(len=88) :: txt ! to store some text
	2060	integer :: ierr,i,j
	2061	integer,parameter :: nblon=180 ! number of longitudes of MONS datasets
	2062	integer,parameter :: nblat=90 ! number of latitudes of MONS datasets
	2063	real :: pi
	2064	real :: longitudes(nblon) ! MONS dataset longitudes
	2065	real :: latitudes(nblat) ! MONS dataset latitudes
	2066	! MONS dataset: mass fraction of H2O where H is assumed to be in H2O
	2067	real :: Hdn(nblon,nblat)
	2068	real :: d21(nblon,nblat)! MONS dataset "depth" (g/cm2)
	2069
	2070	! Extended MONS dataset (for interp_horiz)
	2071	real :: Hdnx(nblon+1,nblat)
	2072	real :: d21x(nblon+1,nblat)
	2073	real :: lon_bound(nblon+1) ! longitude boundaries
	2074	real :: lat_bound(nblat-1) ! latitude boundaries
	2075
	2076	! 1. Initializations:
	2077
	2078	write(,) "Loading MONS data"
	2079
	2080	! Open MONS datafile:
[1918]	2081	open(42,file=trim(datadir)//"/"//trim(filename),
[38]	2082	& status="old",iostat=ierr)
	2083	if (ierr/=0) then
	2084	write(,) "Error in load_MONS_data:"
	2085	write(,) "Failed opening file ",
[1918]	2086	& trim(datadir)//"/"//trim(filename)
	2087	write(,)'1) You can change this directory address in ',
	2088	& 'callfis.def with'
	2089	write(,)' datadir=/path/to/datafiles'
[38]	2090	write(,)'2) If necessary ',trim(filename),
	2091	& ' (and other datafiles)'
	2092	write(,)' can be obtained online at:'
[1383]	2093	write(,)'http://www.lmd.jussieu.fr/~lmdz/planets/mars/datadir'
[38]	2094	CALL ABORT
	2095	else ! skip first line of file (dummy read)
	2096	read(42,*) txt
	2097	endif
	2098
	2099	pi=2.*asin(1.)
	2100
	2101	!2. Load MONS data (on MONS grid)
	2102	do j=1,nblat
	2103	do i=1,nblon
	2104	! swap latitude index so latitudes go from north pole to south pole:
	2105	read(42,*) latitudes(nblat-j+1),longitudes(i),
	2106	& Hdn(i,nblat-j+1),d21(i,nblat-j+1)
	2107	! multiply d21 by 10 to convert from g/cm2 to kg/m2
	2108	d21(i,nblat-j+1)=d21(i,nblat-j+1)*10.0
	2109	enddo
	2110	enddo
	2111	close(42)
	2112
	2113	! there is unfortunately no d21 data for latitudes -77 to -90
	2114	! so we build some by linear interpolation between values at -75
	2115	! and assuming d21=0 at the pole
	2116	do j=84,90 ! latitudes(84)=-77 ; latitudes(83)=-75
	2117	do i=1,nblon
	2118	d21(i,j)=d21(i,83)*((latitudes(j)+90)/15.0)
	2119	enddo
	2120	enddo
	2121
	2122	! 3. Build extended MONS dataset & boundaries (for interp_horiz)
	2123	! longitude boundaries (in radians):
	2124	do i=1,nblon
	2125	! NB: MONS data is every 2 degrees in longitude
	2126	lon_bound(i)=(longitudes(i)+1.0)*pi/180.0
	2127	enddo
	2128	! extra 'modulo' value
	2129	lon_bound(nblon+1)=lon_bound(1)+2.0*pi
	2130
	2131	! latitude boundaries (in radians):
	2132	do j=1,nblat-1
	2133	! NB: Mons data is every 2 degrees in latitude
	2134	lat_bound(j)=(latitudes(j)-1.0)*pi/180.0
	2135	enddo
	2136
	2137	! MONS datasets:
	2138	do j=1,nblat
	2139	Hdnx(1:nblon,j)=Hdn(1:nblon,j)
	2140	Hdnx(nblon+1,j)=Hdnx(1,j)
	2141	d21x(1:nblon,j)=d21(1:nblon,j)
	2142	d21x(nblon+1,j)=d21x(1,j)
	2143	enddo
	2144
	2145	! Interpolate onto GCM grid
	2146	call interp_horiz(Hdnx,MONS_Hdn,nblon,nblat-1,iim,jjm,1,
	2147	& lon_bound,lat_bound,rlonu,rlatv)
	2148	call interp_horiz(d21x,MONS_d21,nblon,nblat-1,iim,jjm,1,
	2149	& lon_bound,lat_bound,rlonu,rlatv)
	2150
	2151	end subroutine

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

Download in other formats:

Original Format