Context Navigation

newstart.F @ 3607

Last change on this file since 3607 was 3539, checked in by tbertrand, 2 months ago
LMDZ.PLUTO: Update and Validation of the Volatile Transport Model (VTM, or "nogcm") Merging missing elements with Pluto.old TB
File size: 148.7 KB

Rev	Line
[3184]	1	C======================================================================
	2	PROGRAM newstart
	3	c=======================================================================
	4	c
	5	c
	6	c Auteur: Christophe Hourdin/Francois Forget/Yann Wanherdrick
	7	c ------
[3371]	8	c Adapted to Pluto: Tanguy Bertrand
	9	c Derniere modif : 06/2024
[3184]	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
	16	use mod_phys_lmdz_para, only: is_parallel, is_sequential,
	17	& is_mpi_root, is_omp_root,
	18	& is_master
	19	use infotrac, only: infotrac_init, nqtot, tname
[3371]	20	USE tracer_h, ONLY: igcm_n2,igcm_ch4_ice,igcm_co_ice,igcm_haze,
	21	& igcm_prec_haze,igcm_co_gas,igcm_ch4_gas,noms
[3184]	22	USE comsoil_h, ONLY: nsoilmx, layer, mlayer, inertiedat
	23	USE surfdat_h, ONLY: phisfi, albedodat,
	24	& zmea, zstd, zsig, zgam, zthe
	25	use datafile_mod, only: datadir, surfdir
	26	use ioipsl_getin_p_mod, only: getin_p
	27	use control_mod, only: day_step, iphysiq, anneeref, planet_type
	28	use phyredem, only: physdem0, physdem1
	29	use iostart, only: open_startphy
	30	use filtreg_mod, only: inifilr
	31	USE mod_const_mpi, ONLY: COMM_LMDZ
	32	USE comvert_mod, ONLY: ap,bp,aps,bps,pa,preff
	33	USE comconst_mod, ONLY: lllm,daysec,dtvr,dtphys,cpp,kappa,
	34	. rad,omeg,g,r,pi
	35	USE serre_mod, ONLY: alphax
	36	USE temps_mod, ONLY: day_ini
	37	USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
	38	use tabfi_mod, only: tabfi
	39	use dimphy, only: init_dimphy
	40	use iniphysiq_mod, only: iniphysiq
	41	use phys_state_var_mod, only: phys_state_var_init
	42	use phyetat0_mod, only: phyetat0
	43	use exner_hyb_m, only: exner_hyb
	44	use geometry_mod, only: longitude, ! longitudes (rad)
	45	& latitude, ! latitudes (rad)
	46	& cell_area ! physics grid area (m2)
	47	implicit none
	48
	49	include "dimensions.h"
	50	integer, parameter :: ngridmx = (2+(jjm-1)*iim - 1/jjm)
	51	include "paramet.h"
	52	include "comgeom2.h"
	53	include "comdissnew.h"
	54	include "netcdf.inc"
	55
	56	c=======================================================================
	57	c Declarations
	58	c=======================================================================
	59
	60	c Variables dimension du fichier "start_archive"
	61	c------------------------------------
	62	CHARACTER relief*3
	63
	64	c Variables pour les lectures NetCDF des fichiers "start_archive"
	65	c--------------------------------------------------
[3382]	66	INTEGER nid_dyn, nid_fi,nid,nvarid,nvarid_input,nvarid_inputs
	67	INTEGER nid_fi_input
[3184]	68	INTEGER length
	69	parameter (length = 100)
	70	INTEGER tab0
	71	INTEGER NB_ETATMAX
	72	parameter (NB_ETATMAX = 100)
	73
	74	REAL date
	75	REAL p_rad,p_omeg,p_g,p_cpp,p_mugaz,p_daysec
	76
	77	c Variable histoire
	78	c------------------
	79	REAL vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm) ! vents covariants
	80	REAL phis(iip1,jjp1)
	81	REAL,ALLOCATABLE :: q(:,:,:,:) ! champs advectes
	82
	83	c autre variables dynamique nouvelle grille
	84	c------------------------------------------
	85	REAL pks(iip1,jjp1)
	86	REAL w(iip1,jjp1,llm+1)
	87	REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
	88	REAL phi(iip1,jjp1,llm)
	89
	90	integer klatdat,klongdat
	91	PARAMETER (klatdat=180,klongdat=360)
	92
	93	c Physique sur grille scalaire
	94	c----------------------------
	95	real zmeaS(iip1,jjp1),zstdS(iip1,jjp1)
	96	real zsigS(iip1,jjp1),zgamS(iip1,jjp1),ztheS(iip1,jjp1)
	97
	98	c variable physique
	99	c------------------
	100	REAL tsurf(ngridmx) ! surface temperature
	101	REAL,ALLOCATABLE :: tsoil(:,:) ! soil temperature
[3539]	102	REAL,ALLOCATABLE :: inertiedat_simple(:,:) ! thermal inertia tmp for dynamico
[3184]	103	REAL emis(ngridmx) ! surface emissivity
	104	real emisread ! added by RW
	105	REAL,ALLOCATABLE :: qsurf(:,:)
[3371]	106	REAL,ALLOCATABLE :: qsurf_tmp(:,:)
[3184]	107	REAL q2(ngridmx,llm+1)
	108	real alb(iip1,jjp1),albfi(ngridmx) ! albedos
	109	real,ALLOCATABLE :: ith(:,:,:),ithfi(:,:) ! thermal inertia (3D)
	110	real surfith(iip1,jjp1),surfithfi(ngridmx) ! surface thermal inertia (2D)
[3371]	111	REAL :: latfi(ngridmx),lonfi(ngridmx),airefi(ngridmx)
[3382]	112	REAL field_input(ngridmx)
	113	REAL,ALLOCATABLE :: field_inputs(:,:)
[3184]	114
[3371]	115	INTEGER i,j,l,isoil,ig,idum,k
[3184]	116	real mugaz ! molar mass of the atmosphere
	117
[3371]	118	integer ierr
[3184]	119
	120	c Variables on the new grid along scalar points
	121	c------------------------------------------------------
	122	REAL t(iip1,jjp1,llm)
	123	REAL tset(iip1,jjp1,llm)
	124	real phisold_newgrid(iip1,jjp1)
[3371]	125	real phisold(iip1,jjp1)
[3184]	126	REAL :: teta(iip1, jjp1, llm)
	127	REAL :: pk(iip1,jjp1,llm)
	128	REAL :: pkf(iip1,jjp1,llm)
	129	REAL :: ps(iip1, jjp1)
	130	REAL :: masse(iip1,jjp1,llm)
	131	REAL :: xpn,xps,xppn(iim),xpps(iim)
	132	REAL :: p3d(iip1, jjp1, llm+1)
	133
	134	c Variable de l'ancienne grille
	135	c------------------------------
	136	real time
	137	real tab_cntrl(100)
	138	real tab_cntrl_bis(100)
	139
	140	c variables diverses
	141	c-------------------
[3371]	142	real choix_1,pp,choice
[3382]	143	integer randtab(ngridmx)
[3184]	144	character*80 fichnom
	145	character*250 filestring
	146	integer Lmodif,iq
	147	character modif*20
	148	real z_reel(iip1,jjp1)
	149	real tsud,albsud,alb_bb,ith_bb,Tiso,Tabove
[3382]	150	real tsurf_bb,tsurf_bb2
	151	real ptoto,pcap,patm,airetot,ptotn,patmn,psea,geop
[3390]	152	real tempsoil(24),levsoil(24)
[3184]	153	character*1 yes
	154	logical :: flagtset=.false. , flagps0=.false.
[3371]	155	real val, val1, val2, val3, val4, dist, ref ! to store temporary variables
	156	real val5, val6, val7, val8, val9, val10,val11, val12 ! to store temporary variables
[3184]	157	real :: iceith=2000 ! thermal inertia of subterranean ice
[3371]	158	integer :: iref,jref,iref1,jref1,iref2,jref2
	159	integer :: igref,igref1,igref2,igref3
[3184]	160
	161	INTEGER :: itau
	162
	163	character(len=20) :: txt ! to store some text
	164	character(len=50) :: surfacefile ! "surface.nc" (or equivalent file)
	165	character(len=150) :: longline
	166	integer :: count
	167	real :: profile(llm+1) ! to store an atmospheric profile + surface value
	168
	169	! added by BC for equilibrium temperature startup
	170	real teque
	171
	172	! added by RW for nuketharsis
	173	real fact1
	174	real fact2
	175
[3371]	176	c Special Pluto Map from Lellouch & Grundy
	177	c------------------------------------------
	178	integer,parameter :: im_plu=360 ! from topo
	179	integer,parameter :: jm_plu=180
	180	real latv_plu(jm_plu+1),lonu_plu(im_plu+1)
	181	real map_pluto_dat(im_plu,jm_plu+1)
[3184]	182
[3371]	183	real N2_pluto_dat(im_plu,jm_plu+1)
	184	real CH4_pluto_dat(im_plu,jm_plu+1)
	185	real CO_pluto_dat(im_plu,jm_plu+1)
	186	real alb_dat(im_plu,jm_plu+1)
	187
	188	real N2_pluto_rein(im_plu+1,jm_plu+1)
	189	real CH4_pluto_rein(im_plu+1,jm_plu+1)
	190	real CO_pluto_rein(im_plu+1,jm_plu+1)
	191	real alb_rein(im_plu+1,jm_plu+1)
	192	real N2_pluto_gcm(iip1,jjp1)
	193	real CH4_pluto_gcm(iip1,jjp1)
	194	real CO_pluto_gcm(iip1,jjp1)
	195	real alb_gcm(iip1,jjp1)
	196
	197	c Special Topo Map mountain
	198	real lat0, lon0
	199	integer ig0
	200	real dist_pluto
	201	real radm,height, phisprev, temp
	202	real preffnew,panew
	203	c Special copy of terrain
	204	real actualmin,angle
	205	integer array_ind(ngridmx)
	206	real array_dist(ngridmx)
	207	real array_angle(ngridmx)
	208
[3184]	209	c sortie visu pour les champs dynamiques
[3198]	210	CALL conf_gcm( 99, .TRUE. )
[3184]	211
	212	cpp = 0.
[3371]	213	preff = 2.
	214	pa = 0.2 ! to ensure disaster rather than minor error if we don`t
[3184]	215	! make deliberate choice of these values elsewhere.
	216
	217	planet_type="generic"
	218
	219	! initialize "serial/parallel" related stuff
	220	! (required because we call tabfi() below, before calling iniphysiq)
	221	is_sequential=.true.
	222	is_parallel=.false.
	223	is_mpi_root=.true.
	224	is_omp_root=.true.
	225	is_master=.true.
	226
	227	! Load tracer number and names:
	228	call infotrac_init
	229	! allocate arrays
	230	allocate(q(iip1,jjp1,llm,nqtot))
	231	allocate(qsurf(ngridmx,nqtot))
[3371]	232	allocate(qsurf_tmp(ngridmx,nqtot))
[3184]	233
	234	! get value of nsoilmx and allocate corresponding arrays
	235	! a default value of nsoilmx is set in comsoil_h
	236	call getin_p("nsoilmx",nsoilmx)
	237
[3539]	238	allocate(inertiedat_simple(ngridmx,nsoilmx))
[3184]	239	allocate(tsoil(ngridmx,nsoilmx))
[3382]	240	allocate(field_inputs(ngridmx,nsoilmx))
[3184]	241	allocate(ith(iip1,jjp1,nsoilmx),ithfi(ngridmx,nsoilmx))
	242
	243	c=======================================================================
	244	c Choice of the start file(s) to use
	245	c=======================================================================
	246	write(,) 'From which kind of files do you want to create new',
	247	. 'start and startfi files'
	248	write(,) ' 0 - from a file start_archive'
	249	write(,) ' 1 - from files start and startfi'
	250
	251	c-----------------------------------------------------------------------
	252	c Open file(s) to modify (start or start_archive)
	253	c-----------------------------------------------------------------------
	254
	255	DO
	256	read(,,iostat=ierr) choix_1
	257	if ((choix_1 /= 0).OR.(choix_1 /=1)) EXIT
	258	ENDDO
	259
	260	c Open start_archive
	261	c ~~~~~~~~~~~~~~~~~~~~~~~~~~
	262	if (choix_1.eq.0) then
	263
	264	write(,) 'Creating start files from:'
	265	write(,) './start_archive.nc'
	266	write(,)
	267	fichnom = 'start_archive.nc'
	268	ierr = NF_OPEN (fichnom, NF_NOWRITE,nid)
	269	IF (ierr.NE.NF_NOERR) THEN
	270	write(6,*)' Problem opening file:',fichnom
	271	write(6,*)' ierr = ', ierr
	272	CALL ABORT
	273	ENDIF
	274	tab0 = 50
	275	Lmodif = 1
	276
	277	c OR open start and startfi files
	278	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	279	else
	280	write(,) 'Creating start files from:'
	281	write(,) './start.nc and ./startfi.nc'
	282	write(,)
	283	fichnom = 'start.nc'
	284	ierr = NF_OPEN (fichnom, NF_NOWRITE,nid_dyn)
	285	IF (ierr.NE.NF_NOERR) THEN
	286	write(6,*)' Problem opening file:',fichnom
	287	write(6,*)' ierr = ', ierr
	288	CALL ABORT
	289	ENDIF
	290
	291	fichnom = 'startfi.nc'
	292	ierr = NF_OPEN (fichnom, NF_NOWRITE,nid_fi)
	293	IF (ierr.NE.NF_NOERR) THEN
	294	write(6,*)' Problem opening file:',fichnom
	295	write(6,*)' ierr = ', ierr
	296	CALL ABORT
	297	ENDIF
	298
	299	tab0 = 0
	300	Lmodif = 0
	301
	302	endif
	303
	304	c=======================================================================
	305	c INITIALISATIONS DIVERSES
	306	c=======================================================================
	307
	308	! Initialize global tracer indexes (stored in tracer.h)
	309	! ... this has to be done before phyetat0
	310	! and requires that "datadir" be correctly initialized
	311	call getin_p("datadir",datadir)
	312	call initracer(ngridmx,nqtot)
	313
	314	! Initialize dimphy module (klon,klev,..)
	315	call init_dimphy(ngridmx,llm)
	316	! Allocate saved arrays (as in firstcall of physiq)
	317	call phys_state_var_init(nqtot)
	318
	319	c-----------------------------------------------------------------------
	320	c Lecture du tableau des parametres du run (pour la dynamique)
	321	c-----------------------------------------------------------------------
	322	if (choix_1.eq.0) then
	323
	324	write(,) 'reading tab_cntrl START_ARCHIVE'
	325	c
	326	ierr = NF_INQ_VARID (nid, "controle", nvarid)
	327	#ifdef NC_DOUBLE
	328	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tab_cntrl)
	329	#else
	330	ierr = NF_GET_VAR_REAL(nid, nvarid, tab_cntrl)
	331	#endif
	332	c
	333	else if (choix_1.eq.1) then
	334
	335	write(,) 'reading tab_cntrl START'
	336	c
	337	ierr = NF_INQ_VARID (nid_dyn, "controle", nvarid)
	338	#ifdef NC_DOUBLE
	339	ierr = NF_GET_VAR_DOUBLE(nid_dyn, nvarid, tab_cntrl)
	340	#else
	341	ierr = NF_GET_VAR_REAL(nid_dyn, nvarid, tab_cntrl)
	342	#endif
	343	c
	344	write(,) 'reading tab_cntrl STARTFI'
	345	c
	346	ierr = NF_INQ_VARID (nid_fi, "controle", nvarid)
	347	#ifdef NC_DOUBLE
	348	ierr = NF_GET_VAR_DOUBLE(nid_fi, nvarid, tab_cntrl_bis)
	349	#else
	350	ierr = NF_GET_VAR_REAL(nid_fi, nvarid, tab_cntrl_bis)
	351	#endif
	352	c
	353	do i=1,50
	354	tab_cntrl(i+50)=tab_cntrl_bis(i)
	355	enddo
	356	write(,) 'printing tab_cntrl', tab_cntrl
	357	do i=1,100
	358	write(,) i,tab_cntrl(i)
	359	enddo
	360
	361	write(,) 'Reading file START'
	362	fichnom = 'start.nc'
	363	CALL dynetat0(fichnom,vcov,ucov,teta,q,masse,
	364	. ps,phis,time)
	365
	366	CALL iniconst
	367	CALL inigeom
	368
	369	! Initialize the physics
	370	CALL iniphysiq(iim,jjm,llm,
	371	& (jjm-1)*iim+2,comm_lmdz,
	372	& daysec,day_ini,dtphys,
	373	& rlatu,rlatv,rlonu,rlonv,
	374	& aire,cu,cv,rad,g,r,cpp,
	375	& 1)
	376
	377	! Lmodif set to 0 to disable modifications possibility in phyeta0
	378	Lmodif=0
	379	write(,) 'Reading file STARTFI'
	380	fichnom = 'startfi.nc'
	381	CALL phyetat0(.true.,ngridmx,llm,fichnom,tab0,Lmodif,nsoilmx,
	382	. nqtot,day_ini,time,
[3539]	383	. tsurf,tsoil,emis,q2,qsurf,inertiedat_simple) !) ! temporary modif by RDW
[3184]	384
	385	! copy albedo and soil thermal inertia on (local) physics grid
	386	do i=1,ngridmx
	387	albfi(i) = albedodat(i)
	388	do j=1,nsoilmx
	389	ithfi(i,j) = inertiedat(i,j)
	390	enddo
	391	! build a surfithfi(:) using 1st layer of ithfi(:), which might
	392	! be needed later on if reinitializing soil thermal inertia
	393	surfithfi(i)=ithfi(i,1)
	394	enddo
	395	! also copy albedo and soil thermal inertia on (local) dynamics grid
	396	! so that options below can manipulate either (but must then ensure
	397	! to correctly recast things on physics grid)
	398	call gr_fi_dyn(1,ngridmx,iip1,jjp1,albfi,alb)
	399	call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
	400	call gr_fi_dyn(1,ngridmx,iip1,jjp1,surfithfi,surfith)
	401
	402	endif
	403	c-----------------------------------------------------------------------
	404	c Initialisation des constantes dynamique
	405	c-----------------------------------------------------------------------
	406
	407	kappa = tab_cntrl(9)
	408	etot0 = tab_cntrl(12)
	409	ptot0 = tab_cntrl(13)
	410	ztot0 = tab_cntrl(14)
	411	stot0 = tab_cntrl(15)
	412	ang0 = tab_cntrl(16)
	413	write(,) "Newstart: kappa,etot0,ptot0,ztot0,stot0,ang0"
	414	write(,) kappa,etot0,ptot0,ztot0,stot0,ang0
	415
	416	! for vertical coordinate
	417	preff=tab_cntrl(18) ! reference surface pressure
	418	pa=tab_cntrl(17) ! reference pressure at which coord is purely pressure
	419	!NB: in start_archive files tab_cntrl(17)=tab_cntrl(18)=0
	420	write(,) "Newstart: preff=",preff," pa=",pa
	421	yes=' '
	422	do while ((yes.ne.'y').and.(yes.ne.'n'))
	423	write(,) "Change the values of preff and pa? (y/n)"
	424	read(*,fmt='(a)') yes
	425	end do
	426	if (yes.eq.'y') then
	427	write(,)"New value of reference surface pressure preff?"
	428	write(,)" (for Mars, typically preff=610)"
	429	read(,) preff
	430	write(,)"New value of reference pressure pa for purely"
	431	write(,)"pressure levels (for hybrid coordinates)?"
	432	write(,)" (for Mars, typically pa=20)"
	433	read(,) pa
	434	endif
	435	c-----------------------------------------------------------------------
	436	c Lecture du tab_cntrl et initialisation des constantes physiques
	437	c - pour start: Lmodif = 0 => pas de modifications possibles
	438	c (modif dans le tabfi de readfi + loin)
	439	c - pour start_archive: Lmodif = 1 => modifications possibles
	440	c-----------------------------------------------------------------------
	441	if (choix_1.eq.0) then
	442	! tabfi requires that input file be first opened by open_startphy(fichnom)
	443	fichnom = 'start_archive.nc'
	444	call open_startphy(fichnom)
	445	call tabfi (ngridmx,nid,Lmodif,tab0,day_ini,lllm,p_rad,
	446	. p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
	447	else if (choix_1.eq.1) then
	448	fichnom = 'startfi.nc'
	449	call open_startphy(fichnom)
	450	Lmodif=1 ! Lmodif set to 1 to allow modifications in phyeta0
	451	call tabfi (ngridmx,nid_fi,Lmodif,tab0,day_ini,lllm,p_rad,
	452	. p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
	453	endif
	454
	455	if (p_omeg .eq. -9999.) then
	456	p_omeg = 8.*atan(1.)/p_daysec
	457	print*,"new value of omega",p_omeg
	458	endif
	459
[3371]	460	c Initialisation coordonnees /aires
	461	c -------------------------------
	462	! Note: rlatu(:) and rlonv(:) are commons defined in "comgeom.h"
	463	! rlatu() and rlonv() are given in radians
	464	latfi(1)=rlatu(1)
	465	lonfi(1)=0.
	466	DO j=2,jjm
	467	DO i=1,iim
	468	latfi((j-2)*iim+1+i)=rlatu(j)
	469	lonfi((j-2)*iim+1+i)=rlonv(i)
	470	ENDDO
	471	ENDDO
	472	latfi(ngridmx)=rlatu(jjp1)
	473	lonfi(ngridmx)=0.
	474
	475	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
	476
[3184]	477	rad = p_rad
	478	omeg = p_omeg
	479	g = p_g
	480	cpp = p_cpp
	481	mugaz = p_mugaz
	482	daysec = p_daysec
	483	kappa = 8.3141000./(p_mugazp_cpp) ! added by RDW
	484
	485	c=======================================================================
	486	c INITIALISATIONS DIVERSES
	487	c=======================================================================
	488	! Load parameters from run.def file
	489	CALL defrun_new( 99, .TRUE. )
	490	! Initialize dynamics geometry and co. (which, when using start.nc may
	491	! have changed because of modifications (tabi, preff, pa) above)
	492	CALL iniconst
	493	CALL inigeom
	494	idum=-1
	495	idum=0
	496
	497	! Initialize the physics for start_archive only
	498	IF (choix_1.eq.0) THEN
	499	CALL iniphysiq(iim,jjm,llm,
	500	& (jjm-1)*iim+2,comm_lmdz,
	501	& daysec,day_ini,dtphys,
	502	& rlatu,rlatv,rlonu,rlonv,
	503	& aire,cu,cv,rad,g,r,cpp,
	504	& 1)
	505	ENDIF
	506
	507	c=======================================================================
	508	c lecture topographie, albedo, inertie thermique, relief sous-maille
	509	c=======================================================================
	510
	511	if (choix_1.eq.0) then ! for start_archive files,
	512	! where an external "surface.nc" file is needed
	513
	514	c do while((relief(1:3).ne.'mol').AND.(relief(1:3).ne.'pla'))
	515	c write(,)
	516	c write(,) 'choix du relief (mola,pla)'
	517	c write(,) '(Topographie MGS MOLA, plat)'
	518	c read(*,fmt='(a3)') relief
	519	relief="mola"
	520	c enddo
	521
	522	! First get the correct value of datadir, if not already done:
	523	! default 'datadir' is set in "datafile_mod"
	524	call getin_p("datadir",datadir)
	525	write(,) 'Available surface data files are:'
	526	filestring='ls '//trim(datadir)//'/'//
	527	& trim(surfdir)//' \| grep .nc'
	528	call system(filestring)
	529	! but in ye old days these files were in datadir, so scan it as well
	530	! for the sake of retro-compatibility
	531	filestring='ls '//trim(datadir)//'/'//' \| grep .nc'
	532	call system(filestring)
	533
	534	write(,) ''
	535	write(,) 'Please choose the relevant file',
	536	& ' (e.g. "surface_mars.nc")'
	537	write(,) ' or "none" to not use any (and set planetary'
	538	write(,) ' albedo and surface thermal inertia)'
	539	read(*,fmt='(a50)') surfacefile
	540
	541	if (surfacefile.ne."none") then
	542
	543	CALL datareadnc(relief,surfacefile,phis,alb,surfith,
	544	& zmeaS,zstdS,zsigS,zgamS,ztheS)
	545	else
	546	! specific case when not using a "surface.nc" file
	547	phis(:,:)=0
	548	zmeaS(:,:)=0
	549	zstdS(:,:)=0
	550	zsigS(:,:)=0
	551	zgamS(:,:)=0
	552	ztheS(:,:)=0
	553
	554	write(,) "Enter value of albedo of the bare ground:"
	555	read(,) alb(1,1)
	556	alb(:,:)=alb(1,1)
	557
	558	write(,) "Enter value of thermal inertia of soil:"
	559	read(,) surfith(1,1)
	560	surfith(:,:)=surfith(1,1)
	561
	562	endif ! of if (surfacefile.ne."none")
	563
	564	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
	565	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,surfith,surfithfi)
	566	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
	567	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zmeaS,zmea)
	568	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zstdS,zstd)
	569	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zsigS,zsig)
	570	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zgamS,zgam)
	571	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ztheS,zthe)
	572
	573	endif ! of if (choix_1.eq.0)
	574
	575	c=======================================================================
	576	c Lecture des fichiers (start ou start_archive)
	577	c=======================================================================
	578
	579	if (choix_1.eq.0) then
	580
	581	write(,) 'Reading file START_ARCHIVE'
	582	CALL lect_start_archive(ngridmx,llm,
	583	& date,tsurf,tsoil,emis,q2,
	584	& t,ucov,vcov,ps,teta,phisold_newgrid,q,qsurf,
[3198]	585	& surfith,nid)
	586
[3184]	587	write(,) "OK, read start_archive file"
	588	! copy soil thermal inertia
	589	ithfi(:,:)=inertiedat(:,:)
	590
	591	ierr= NF_CLOSE(nid)
	592
	593	else if (choix_1.eq.1) then
	594	!do nothing, start and startfi have already been read
	595	else
	596	CALL exit(1)
	597	endif
	598
	599	dtvr = daysec/FLOAT(day_step)
	600	dtphys = dtvr * FLOAT(iphysiq)
	601
	602	c=======================================================================
	603	c
	604	c=======================================================================
	605
	606	do ! infinite loop on list of changes
	607
[3371]	608	write(,)
	609	write(,)
	610	write(,) 'List of possible changes :'
	611	write(,) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
	612	write(,)
	613	write(,) 'flat : no topography ("aquaplanet")'
	614	write(,) 'set_ps_to_preff : used if changing preff with topo'
	615	write(,) 'qname : change tracer name'
	616	write(,) 't=profile : read temperature profile in profile.in'
	617	write(,) 'q=0 : ALL tracer =zero'
	618	write(,) 'q=x : give a specific uniform value to one tracer'
	619	write(,) 'qs=x : give a uniform value to a surface tracer'
	620	write(,) 'q=profile : specify a profile for a tracer'
	621	write(,) 'qnogcm: initial tracer for GCM from nogcm (CH4,CO)'
	622	write(,) 'isotherm : Isothermal Temperatures'
	623	write(,) 'tempprof : specific Temperature profile'
	624	write(,) 'initsoil : initialize soil Temperatures'
	625	write(,) 'ptot : change total pressure'
	626	write(,) 'therm_ini_s: set soil TI to reference surf. values'
	627	write(,) 'inert3d: give uniform value of thermal inertia'
	628	write(,) 'subsoilice_n: put deep underground ice in N. hemis'
	629	write(,) 'subsoilice_s: put deep underground ice in S. hemis'
	630	write(,) 'reservoir: increase/decrease reservoir of ice'
	631	write(,) 'reservoir_SP: increase/decrease reservoir in SP'
	632	write(,) 'plutomap: initialize surface like pluto from ...'
	633	write(,) 'subsoil_n2: diu-sea thermal inertia for N2 only'
	634	write(,) 'subsoil_ch4: diu-sea thermal inertia for CH4 only'
	635	write(,) 'subsoil_all: diu-sea thermal inertia for all terr'
	636	write(,) 'subsoil_alb: diu-sea thermal inertia from albedo'
	637	write(,) 'diurnal_TI: diurnal thermal inertia for all terr'
	638	write(,) 'initsurf: initial surface temperature (global)'
	639	write(,) 'modtsurf: initial surface temperature (global)'
	640	write(,) 'copylat: copy tsurf and tsoil latitude'
	641	write(,) 'copylatlon: copy tsurf and tsoil lat/lon'
	642	write(,) 'copylon: copy tsurf tsoil lat, n2surf and phisfi'
	643	write(,) 'tsurfice: temperature depending on N2 ice'
	644	write(,) 'icarus: option to change surface/soil temperature'
	645	write(,) 'icarus_ch4: special option to add ch4'
	646	write(,) 'delfrostch4: delete ch4 frost'
	647	write(,) 'modch4: remove/modify amount ch4 frost'
	648	write(,) 'modch4n2: modify amount ch4 frost according N2'
	649	write(,) 'modco: remove/modify amount co frost'
	650	write(,) 'modn2: remove/modify amount n2 ice'
	651	write(,) 'modcoch4: remove/modify co ch4 where no n2 '
	652	write(,) 'checktsoil: change tsoil where no n2 '
	653	write(,) 'non2noco: if no n2, no co '
	654	write(,) 'modn2_topo: modify n2 ice topo according to topo'
	655	write(,) 'modwheren2: modify n2 where already n2'
	656	write(,) 'modn2HD: modify n2 for HD runs'
	657	write(,) 'modn2HD_SP: modify n2 for HD runs in SP'
	658	write(,) 'globch4co: add/remove global amount ch4co frost'
	659	write(,) 'source_ch4: add source ch4'
	660	write(,) 'nomountain: remove pluto mountains (!)'
	661	write(,) 'relief: add pluto crater or mountain'
	662	write(,) 'seticeNH: set ice for initial start with NH topo'
	663	write(,) 'topo_sp: change topography of SP'
	664	write(,) 'fill_sp: fill sp with N2 ice and adjust phisfi'
	665	write(,) 'fill_sp_inv: fill inverted sp and adjust'
	666	write(,) 'subsoil_SP: diu-sea thermal inertia for SP '
	667	write(,) 'bladed: add ch4 on bladed terrains'
	668	write(,) 'cpbladed: add extension bladed terrains'
	669	write(,) 'slope: add slope over all long. (for triton)'
	670	write(,) 'digsp: specific to dig SP'
	671	write(,) 'bugn2: correct bug of warm n2 due to HighRes'
	672	write(,) 'bugch4: correct bug of warm ch4 due to HighRes'
	673	write(,) 'flatnosp : topo flat except SP (topo controled)'
	674	write(,) 'flatregion: topo flat for specific region'
	675	write(,) 'changetopo: change topo'
	676	write(,) 'asymtopo: N-S asym topo in tanh'
	677	write(,) 'corrtopo: correction topo bug'
	678	write(,) 'adjustphi: adjust phisfi according to N2 mass'
	679	write(,) 'correctphi: adjust phisfi'
	680	write(,) 'correctps: test to correct ps'
	681	write(,) 'toponoise: no flat topo'
	682	write(,) 'asymtriton: asymetry in longitude for triton'
	683	write(,) 'copytsoil: copy 2D tsoil field from nc file'
	684	write(,) 'albedomap: read in an albedomap albedo.nc'
	685	write(,) 'mod_distrib_ice: modify ice distrib. from albedo'
	686
	687	write(,)
	688	write(,) 'Please note that emis and albedo are set in physiq'
	689	write(,)
[3184]	690
	691	write(,)
	692	write(,) 'Change to perform ?'
	693	write(,) ' (enter keyword or return to end)'
	694	write(,)
	695
	696	read(*,fmt='(a20)') modif
	697	if (modif(1:1) .eq. ' ')then
	698	exit ! exit loop on changes
	699	endif
	700
	701	write(,)
	702	write(,) trim(modif) , ' : '
	703
	704	c 'flat : no topography ("aquaplanet")'
	705	c -------------------------------------
	706	if (trim(modif) .eq. 'flat') then
	707	c set topo to zero
	708	z_reel(1:iip1,1:jjp1)=0
	709	phis(1:iip1,1:jjp1)=g*z_reel(1:iip1,1:jjp1)
	710	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
	711	write(,) 'topography set to zero.'
	712	write(,) 'WARNING : the subgrid topography parameters',
	713	& ' were not set to zero ! => set calllott to F'
	714
	715	c Choice of surface pressure
	716	yes=' '
	717	do while ((yes.ne.'y').and.(yes.ne.'n'))
	718	write(,) 'Do you wish to choose homogeneous surface',
	719	& 'pressure (y) or let newstart interpolate ',
	720	& ' the previous field (n)?'
	721	read(*,fmt='(a)') yes
	722	end do
	723	if (yes.eq.'y') then
	724	flagps0=.true.
	725	write(,) 'New value for ps (Pa) ?'
	726	201 read(,,iostat=ierr) patm
	727	if(ierr.ne.0) goto 201
	728	write(,) patm
	729	if (patm.eq.-9999.) then
	730	patm = preff
	731	write(,) "we set patm = preff", preff
	732	endif
	733	write(,)
	734	write(,) ' new ps everywhere (Pa) = ', patm
	735	write(,)
	736	do j=1,jjp1
	737	do i=1,iip1
	738	ps(i,j)=patm
	739	enddo
	740	enddo
	741	end if
	742
	743	c 'set_ps_to_preff' : used if changing preff with topo
	744	c ----------------------------------------------------
	745	else if (trim(modif) .eq. 'set_ps_to_preff') then
	746	do j=1,jjp1
	747	do i=1,iip1
	748	ps(i,j)=preff
	749	enddo
	750	enddo
	751
[3382]	752	c ptot : Modification of the total pressure: ice + current atmosphere
	753	c -------------------------------------------------------------------
	754	else if (modif(1:len_trim(modif)).eq.'ptot') then
	755
	756	c calcul de la pression totale glace + atm actuelle
	757	patm=0.
	758	airetot=0.
	759	pcap=0.
	760	DO j=1,jjp1
	761	DO i=1,iim
	762	ig=1+(j-2)*iim +i
	763	if(j.eq.1) ig=1
	764	if(j.eq.jjp1) ig=ngridmx
	765	patm = patm + ps(i,j)*aire(i,j)
	766	airetot= airetot + aire(i,j)
	767	ENDDO
	768	ENDDO
	769	ptoto = pcap + patm
	770
	771	print*,'Current total pressure at surface ',
	772	& ptoto/airetot
	773
	774	print*,'new value?'
	775	read(,) ptotn
	776	ptotn=ptotn*airetot
	777	patmn=ptotn-pcap
	778	print*,'ptoto,patm,ptotn,patmn'
	779	print*,ptoto,patm,ptotn,patmn
	780	print*,'Mult. factor for pressure (atm only)', patmn/patm
	781	do j=1,jjp1
	782	do i=1,iip1
	783	ps(i,j)=ps(i,j)*patmn/patm
	784	enddo
	785	enddo
	786
	787	c Correction pour la conservation des traceurs
	788	yes=' '
	789	do while ((yes.ne.'y').and.(yes.ne.'n'))
	790	write(,) 'Do you wish to conserve tracer total mass (y)',
	791	& ' or tracer mixing ratio (n) ?'
	792	read(*,fmt='(a)') yes
	793	end do
	794
	795	if (yes.eq.'y') then
	796	write(,) 'OK : conservation of tracer total mass'
	797	DO iq =1, nqtot
	798	DO l=1,llm
	799	DO j=1,jjp1
	800	DO i=1,iip1
	801	q(i,j,l,iq)=q(i,j,l,iq)*patm/patmn
	802	ENDDO
	803	ENDDO
	804	ENDDO
	805	ENDDO
	806	else
	807	write(,) 'OK : conservation of tracer mixing ratio'
	808	end if
	809
[3184]	810	c qname : change tracer name
	811	c --------------------------
	812	else if (trim(modif).eq.'qname') then
	813	yes='y'
	814	do while (yes.eq.'y')
	815	write(,) 'Which tracer name do you want to change ?'
	816	do iq=1,nqtot
	817	write(*,'(i3,a3,a20)')iq,' : ',trim(tname(iq))
	818	enddo
	819	write(*,'(a35,i3)')
	820	& '(enter tracer number; between 1 and ',nqtot
	821	write(,)' or any other value to quit this option)'
	822	read(,) iq
	823	if ((iq.ge.1).and.(iq.le.nqtot)) then
	824	write(,)'Change tracer name ',trim(tname(iq)),' to ?'
	825	read(,) txt
	826	tname(iq)=txt
	827	write(,)'Do you want to change another tracer name (y/n)?'
	828	read(*,'(a)') yes
	829	else
	830	! inapropiate value of iq; quit this option
	831	yes='n'
	832	endif ! of if ((iq.ge.1).and.(iq.le.nqtot))
	833	enddo ! of do while (yes.ne.'y')
	834
	835	c q=0 : set tracers to zero
	836	c -------------------------
	837	else if (trim(modif).eq.'q=0') then
	838	c mise a 0 des q (traceurs)
	839	write(,) 'Tracers set to 0 (1.E-30 in fact)'
	840	DO iq =1, nqtot
	841	DO l=1,llm
	842	DO j=1,jjp1
	843	DO i=1,iip1
	844	q(i,j,l,iq)=1.e-30
	845	ENDDO
	846	ENDDO
	847	ENDDO
	848	ENDDO
	849
	850	c set surface tracers to zero
	851	DO iq =1, nqtot
	852	DO ig=1,ngridmx
	853	qsurf(ig,iq)=0.
	854	ENDDO
	855	ENDDO
	856
	857	c q=x : initialise tracer manually
	858	c --------------------------------
	859	else if (trim(modif).eq.'q=x') then
	860	write(,) 'Which tracer do you want to modify ?'
	861	do iq=1,nqtot
	862	write(,)iq,' : ',trim(tname(iq))
	863	enddo
	864	write(,) '(choose between 1 and ',nqtot,')'
	865	read(,) iq
	866	write(,)'mixing ratio of tracer ',trim(tname(iq)),
	867	& ' ? (kg/kg)'
	868	read(,) val
	869	DO l=1,llm
	870	DO j=1,jjp1
	871	DO i=1,iip1
	872	q(i,j,l,iq)=val
	873	ENDDO
	874	ENDDO
	875	ENDDO
[3390]	876	c write(,) 'SURFACE value of tracer ',trim(tname(iq)),
	877	c & ' ? (kg/m2)'
	878	c read(,) val
	879	c DO ig=1,ngridmx
	880	c qsurf(ig,iq)=val
	881	c ENDDO
[3184]	882
	883	c qs=x : initialise surface tracer manually
	884	c --------------------------------
	885	else if (trim(modif).eq.'qs=x') then
	886	write(,) 'Which tracer do you want to modify ?'
	887	do iq=1,nqtot
	888	write(,)iq,' : ',trim(tname(iq))
	889	enddo
	890	write(,) '(choose between 1 and ',nqtot,')'
	891	read(,) iq
	892	write(,) 'SURFACE value of tracer ',trim(tname(iq)),
	893	& ' ? (kg/m2)'
	894	read(,) val
	895	DO ig=1,ngridmx
	896	qsurf(ig,iq)=val
	897	ENDDO
	898
	899	c t=profile : initialize temperature with a given profile
	900	c -------------------------------------------------------
	901	else if (trim(modif) .eq. 't=profile') then
	902	write(,) 'Temperature profile from ASCII file'
	903	write(,) "'profile.in' e.g. 1D output"
	904	write(,) "(one value per line in file; starting with"
	905	write(,) "surface value, the 1st atmospheric layer"
	906	write(,) "followed by 2nd, etc. up to top of atmosphere)"
	907	txt="profile.in"
	908	open(33,file=trim(txt),status='old',form='formatted',
	909	& iostat=ierr)
	910	if (ierr.eq.0) then
	911	! OK, found file 'profile_...', load the profile
	912	do l=1,llm+1
	913	read(33,*,iostat=ierr) profile(l)
	914	write(,) profile(l)
	915	if (ierr.ne.0) then ! something went wrong
	916	exit ! quit loop
	917	endif
	918	enddo
	919	if (ierr.eq.0) then
	920	tsurf(1:ngridmx)=profile(1)
	921	tsoil(1:ngridmx,1:nsoilmx)=profile(1)
	922	do l=1,llm
	923	Tset(1:iip1,1:jjp1,l)=profile(l+1)
	924	flagtset=.true.
	925	enddo
	926	ucov(1:iip1,1:jjp1,1:llm)=0.
	927	vcov(1:iip1,1:jjm,1:llm)=0.
	928	q2(1:ngridmx,1:llm+1)=0.
	929	else
	930	write(,)'problem reading file ',trim(txt),' !'
	931	write(,)'No modifications to temperature'
	932	endif
	933	else
	934	write(,)'Could not find file ',trim(txt),' !'
	935	write(,)'No modifications to temperature'
	936	endif
	937
	938	c q=profile : initialize tracer with a given profile
	939	c --------------------------------------------------
	940	else if (trim(modif) .eq. 'q=profile') then
	941	write(,) 'Tracer profile will be sought in ASCII file'
	942	write(,) "'profile_tracer' where 'tracer' is tracer name"
	943	write(,) "(one value per line in file; starting with"
	944	write(,) "surface value, the 1st atmospheric layer"
	945	write(,) "followed by 2nd, etc. up to top of atmosphere)"
	946	write(,) 'Which tracer do you want to set?'
	947	do iq=1,nqtot
	948	write(,)iq,' : ',trim(tname(iq))
	949	enddo
	950	write(,) '(choose between 1 and ',nqtot,')'
	951	read(,) iq
	952	if ((iq.lt.1).or.(iq.gt.nqtot)) then
	953	! wrong value for iq, go back to menu
	954	write(,) "wrong input value:",iq
	955	cycle
	956	endif
	957	! look for input file 'profile_tracer'
	958	txt="profile_"//trim(tname(iq))
	959	open(41,file=trim(txt),status='old',form='formatted',
	960	& iostat=ierr)
	961	if (ierr.eq.0) then
	962	! OK, found file 'profile_...', load the profile
	963	do l=1,llm+1
	964	read(41,*,iostat=ierr) profile(l)
	965	if (ierr.ne.0) then ! something went wrong
	966	exit ! quit loop
	967	endif
	968	enddo
	969	if (ierr.eq.0) then
	970	! initialize tracer values
	971	qsurf(:,iq)=profile(1)
	972	do l=1,llm
	973	q(:,:,l,iq)=profile(l+1)
	974	enddo
	975	write(,)'OK, tracer ',trim(tname(iq)),' initialized '
	976	& ,'using values from file ',trim(txt)
	977	else
	978	write(,)'problem reading file ',trim(txt),' !'
	979	write(,)'No modifications to tracer ',trim(tname(iq))
	980	endif
	981	else
	982	write(,)'Could not find file ',trim(txt),' !'
	983	write(,)'No modifications to tracer ',trim(tname(iq))
	984	endif
	985
[3371]	986	c qnogcm : initialise tracer from nogcm (CH4, CO)
	987	c --------------------------------
	988	else if (modif(1:len_trim(modif)).eq.'qnogcm') then
	989	write(,) 'Which tracer do you want to modify ?'
	990	write(,) 'Should be ch4_gas and co_gas'
	991	do iq=1,nqtot
	992	write(,)iq,' : ',trim(noms(iq)),' : ',q(1,1,1,iq)
	993	enddo
	994	write(,) '(choose between 1 and ',nqtot,')'
	995	read(,) iq
	996	DO l=1,llm
	997	DO j=1,jjp1
	998	DO i=1,iip1
	999	q(i,j,l,iq)=q(i,j,1,iq)
	1000	ENDDO
	1001	ENDDO
	1002	ENDDO
[3198]	1003
[3371]	1004	c isothermal temperatures
	1005	c ------------------------------------------------
	1006	else if (modif(1:len_trim(modif)) .eq. 'isotherm') then
[3228]	1007
[3371]	1008	write(,)'Isothermal temperature of the atmosphere'
	1009	write(,) 'Value of THIS temperature ? :'
	1010	203 read(,,iostat=ierr) Tset(1,1,1)
	1011	if(ierr.ne.0) goto 203
	1012	flagtset=.true.
	1013	DO l=1,llm
	1014	DO j=1,jjp1
	1015	DO i=1,iip1
	1016	Tset(i,j,l)=Tset(1,1,1)
	1017	ENDDO
	1018	ENDDO
	1019	ENDDO
	1020	write(,) 'atmospheric temp =', Tset(2,2,2)
[3228]	1021
[3371]	1022	c specific temperature profile
	1023	c ------------------------------------------------
	1024	else if (modif(1:len_trim(modif)) .eq. 'tempprof') then
[3228]	1025
[3371]	1026	write(,) 'phi='
	1027	write(,) phi(1,1,:)
	1028	write(,) 'temperature profile in the atmosphere'
	1029	write(,) 'Value of THIS temperature ? :'
	1030	204 read(,,iostat=ierr) Tset(1,1,1)
	1031	if(ierr.ne.0) goto 204
	1032	flagtset=.true.
	1033	DO l=1,llm
	1034	DO j=1,jjp1
	1035	DO i=1,iip1
	1036	Tset(i,j,l)=Tset(1,1,1)
	1037	ENDDO
	1038	ENDDO
	1039	ENDDO
	1040	write(,) 'atmospheric temp =', Tset(2,2,2)
[3228]	1041
[3371]	1042	c initsoil: subsurface temperature
	1043	c ---------------------------------
	1044	else if (modif(1:len_trim(modif)) .eq. 'initsoil') then
[3228]	1045
[3371]	1046	write(,)'Isothermal temperature of the subsurface'
	1047	write(,) 'Value of this temperature ? :'
	1048	303 read(,,iostat=ierr) Tiso
	1049	if(ierr.ne.0) goto 303
	1050
	1051	do l=1,nsoilmx
	1052	do ig=1, ngridmx
	1053	tsoil(ig,l) = Tiso
	1054	end do
	1055	end do
	1056
	1057	c icarus : changing tsoil tsurf on latitudinal bands
	1058	c -----------------------------------------------------
	1059	else if (modif(1:len_trim(modif)) .eq. 'icarus') then
	1060
	1061	write(,) 'At which lat lon do you want to extract the
	1062	& reference tsurf/tsoil profile ?'
	1063	407 read(,,iostat=ierr) val,val2
	1064	if(ierr.ne.0) goto 407
	1065	write(,) 'You want lat =',val
	1066	write(,) 'You want lon =',val2
	1067	dist=0
	1068	ref=1000
	1069	do j=1,ngridmx-1
	1070	dist=sqrt((lonfi(j)180./pi-val2)*2+
	1071	& (latfi(j)180./pi-val)*2)
	1072	IF (dist.lt.ref) then
	1073	ref=dist
	1074	jref=j
	1075	ENDIF
	1076	ENDDO
	1077
	1078	write(,)'Will use nearest grid point which is:',
	1079	& latfi(jref)180./pi,lonfi(jref)180./pi
	1080
	1081	! Extraction of the profile
	1082	write(,) 'Profile is =',tsoil(jref,:)
	1083	write(,) 'tsurf is =',tsurf(jref)
	1084	write(,) 'Choice lat for latitudinal band with this profile'
	1085	408 read(,,iostat=ierr) val3
	1086	if(ierr.ne.0) goto 408
	1087	write(,) 'Choice delta lat for this latitudinal band'
	1088	409 read(,,iostat=ierr) val4
	1089	if(ierr.ne.0) goto 409
	1090	write(,) 'Choice delta temp (K) for this latitudinal band'
	1091	410 read(,,iostat=ierr) val5
	1092	if(ierr.ne.0) goto 410
	1093	write(,) 'How much N2 ice should I put on this band (kg/m2)'
	1094	415 read(,,iostat=ierr) choice
	1095	if(ierr.ne.0) goto 415
	1096	write(,) 'Choice lat2'
	1097	411 read(,,iostat=ierr) val6
	1098	if(ierr.ne.0) goto 411
	1099	write(,) 'Choice delta lat for this latitudinal band'
	1100	412 read(,,iostat=ierr) val7
	1101	if(ierr.ne.0) goto 412
	1102	write(,) 'Choice delta temp (K) for this latitudinal band'
	1103	413 read(,,iostat=ierr) val8
	1104	if(ierr.ne.0) goto 413
	1105
[3228]	1106	DO ig=1,ngridmx
[3371]	1107	IF ( (latfi(ig)*180./pi.ge.(val3-val4/2.)) .and.
	1108	& (latfi(ig)*180./pi.le.(val3+val4/2.)) .and.
	1109	& (qsurf(ig,igcm_n2).lt.0.001) ) then
	1110	tsurf(ig)=tsurf(jref)+val5
	1111	DO l=1,nsoilmx
	1112	tsoil(ig,l)=tsoil(jref,l)+val5
	1113	ENDDO
	1114	qsurf(ig,igcm_n2)=choice
	1115	ENDIF
	1116	IF ( (latfi(ig)*180./pi.ge.(val6-val7/2.)) .and.
	1117	& (latfi(ig)*180./pi.le.(val6+val7/2.)) .and.
	1118	& (qsurf(ig,igcm_n2).lt.0.001) ) then
	1119	tsurf(ig)=tsurf(jref)+val8
	1120	DO l=1,nsoilmx
	1121	tsoil(ig,l)=tsoil(jref,l)+val8
	1122	ENDDO
	1123	ENDIF
	1124	ENDDO
	1125
	1126	c icarus_ch4 : changing tsoil tsurf and adding ch4
	1127	c -----------------------------------------------------
	1128	else if (modif(1:len_trim(modif)) .eq. 'icarus_ch4') then
	1129
	1130	write(,) 'At which lat lon do you want to extract the
	1131	& reference tsurf/tsoil profile ?'
	1132	416 read(,,iostat=ierr) val,val2
	1133	if(ierr.ne.0) goto 416
	1134	write(,) 'You want lat =',val
	1135	write(,) 'You want lon =',val2
	1136	dist=0
	1137	ref=1000
	1138	do j=1,ngridmx-1
	1139	dist=sqrt((lonfi(j)180./pi-val2)*2+
	1140	& (latfi(j)180./pi-val)*2)
	1141	IF (dist.lt.ref) then
	1142	ref=dist
	1143	jref=j
	1144	ENDIF
	1145	ENDDO
	1146
	1147	write(,)'Will use nearest grid point which is:',
	1148	& latfi(jref)180./pi,lonfi(jref)180./pi
	1149
	1150	! Extraction of the profile
	1151	write(,) 'Profile is =',tsoil(jref,:)
	1152	write(,) 'tsurf is =',tsurf(jref)
	1153	write(,) 'Choice band : lat1 and lat2 ?'
	1154	417 read(,,iostat=ierr) val3,val4
	1155	if(ierr.ne.0) goto 417
	1156	write(,) 'Choice temp coefficient for this latitudinal band'
	1157	418 read(,,iostat=ierr) val5
	1158	if(ierr.ne.0) goto 418
	1159	write(,) 'How much CH4 ice do I put on this band (kg/m2)'
	1160	419 read(,,iostat=ierr) choice
	1161	if(ierr.ne.0) goto 419
	1162
	1163	DO ig=1,ngridmx
	1164	IF ( (latfi(ig)*180./pi.ge.val3) .and.
	1165	& (latfi(ig)*180./pi.le.val4) .and.
	1166	& (qsurf(ig,igcm_ch4_ice).lt.0.001) ) then
	1167	tsurf(ig)=tsurf(jref)*val5
	1168	DO l=1,nsoilmx
	1169	tsoil(ig,l)=tsoil(jref,l)*val5
	1170	ENDDO
	1171	qsurf(ig,igcm_ch4_ice)=choice
	1172	ENDIF
	1173	ENDDO
	1174
	1175	c globch4co : adding/remove global ch4 co
	1176	c ----------------------------------------------------------
	1177	else if (modif(1:len_trim(modif)) .eq. 'globch4co') then
	1178
	1179	write(,) 'Adding or removing ch4 co '
	1180	write(,) 'Choice amount add/less ch4 ice (0 = rm all)'
	1181	431 read(,,iostat=ierr) val3
	1182	if(ierr.ne.0) goto 431
	1183	write(,) 'Choice amount add/less co ice (0 = rm all)'
	1184	432 read(,,iostat=ierr) val4
	1185	if(ierr.ne.0) goto 432
	1186	IF (val3.eq.0.) then
	1187	DO ig=1,ngridmx
	1188	qsurf(ig,igcm_ch4_ice)=0.
[3228]	1189	ENDDO
[3371]	1190	ELSE
	1191	DO ig=1,ngridmx
	1192	qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)+val3
	1193	ENDDO
	1194	ENDIF
	1195	IF (val4.eq.0.) then
	1196	DO ig=1,ngridmx
	1197	qsurf(ig,igcm_co_ice)=0.
	1198	ENDDO
	1199	ELSE
	1200	DO ig=1,ngridmx
	1201	qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)+val4
	1202	ENDDO
	1203	ENDIF
	1204
	1205	c bladed : adding/remove ch4 on bladed terrains
	1206	c ----------------------------------------------------------
	1207	else if (modif(1:len_trim(modif)).eq.'bladed') then
	1208
	1209	write(,) 'Adding or removing ch4 on the bladed terrains'
	1210	write(,) 'Choice band : lat1 and lat2 ?'
	1211	450 read(,,iostat=ierr) val,val2
	1212	if(ierr.ne.0) goto 450
	1213	write(,) 'Choice band : lon1 and lon2 ?'
	1214	451 read(,,iostat=ierr) val3,val4
	1215	if(ierr.ne.0) goto 451
	1216	write(,) 'Choice phisfi minimum ?'
	1217	454 read(,,iostat=ierr) choice
	1218	if(ierr.ne.0) goto 454
	1219	write(,) 'Choice multiplicative factor'
	1220	452 read(,,iostat=ierr) val5
	1221	if(ierr.ne.0) goto 452
	1222	write(,) 'Choice additional ch4'
	1223	453 read(,,iostat=ierr) val6
	1224	if(ierr.ne.0) goto 453
	1225	write(,) 'Choice index for tsurf tsoil'
	1226	449 read(,,iostat=ierr) iref
	1227	if(ierr.ne.0) goto 449
	1228
	1229	! shift
	1230	DO ig=1,ngridmx
	1231	IF ( (latfi(ig)*180./pi.ge.val) .and.
	1232	& (latfi(ig)*180./pi.le.val2) .and.
	1233	& (lonfi(ig)*180./pi.ge.val3) .and.
	1234	& (lonfi(ig)*180./pi.le.val4) .and.
	1235	& (phisfi(ig).gt.choice) ) then
	1236	qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)*val5
	1237	qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)+val6
	1238	tsurf(ig)=tsurf(iref)
	1239	DO l=1,nsoilmx
	1240	tsoil(ig,l)=tsoil(iref,l)
	1241	ENDDO
	1242	ENDIF
[3228]	1243	ENDDO
	1244
[3371]	1245	c cpbladed : add extension bladed terrains
	1246	c ----------------------------------------------------------
	1247	else if (modif(1:len_trim(modif)).eq.'cpbladed') then
[3228]	1248
[3371]	1249	write(,) 'Choice lat,lon, centre of band to be copied ?'
	1250	560 read(,,iostat=ierr) val,val2
	1251	if(ierr.ne.0) goto 560
	1252	write(,) 'Choice distance (km) from there defining the area'
	1253	561 read(,,iostat=ierr) val3
	1254	if(ierr.ne.0) goto 561
	1255	write(,) 'Nb of copies ?'
	1256	562 read(,,iostat=ierr) l
	1257	if(ierr.ne.0) goto 562
[3228]	1258
[3371]	1259	! Get index correponding to central points
	1260	! 1/ Reference point
	1261	igref=1
	1262	actualmin=1000.
	1263	do ig=1,ngridmx
	1264	val6=dist_pluto(latfi(ig),lonfi(ig),
	1265	& valpi/180.,val2pi/180.)
	1266	if (val6.lt.actualmin) then
	1267	actualmin=val6
	1268	igref=ig
	1269	endif
	1270	enddo
[3228]	1271
[3371]	1272	do k=1,l
	1273
	1274	write(,) 'Choice lat,lon where terrain is copied'
	1275	563 read(,,iostat=ierr) val4,val5
	1276	if(ierr.ne.0) goto 563
	1277
	1278	if (val5.gt.180.) then
	1279	val5=val5-360.
	1280	endif
	1281
	1282	! 2/ Target point
	1283	igref2=1
	1284	actualmin=1000.
	1285	do ig=1,ngridmx
	1286	val6=dist_pluto(latfi(ig),lonfi(ig),
	1287	& val4pi/180.,val5pi/180.)
	1288	if (val6.lt.actualmin) then
	1289	actualmin=val6
	1290	igref2=ig
	1291	endif
	1292	enddo
	1293
	1294	! for each point within A1, get distance and angle
	1295	! save in a table
	1296	i=1
	1297	array_ind(:)=0
	1298	array_dist(:)=1000.
	1299	array_angle(:)=0.
	1300	do ig=1,ngridmx
	1301	val6=dist_pluto(latfi(ig),lonfi(ig),
	1302	& valpi/180.,val2pi/180.)
	1303	if (val6.lt.val3) then
	1304	array_ind(i)=ig
	1305	array_dist(i)=val6
	1306	array_angle(i)=atan2(val-latfi(ig)*180./pi,
	1307	& val2-lonfi(ig)*180./pi)
	1308	i=i+1
	1309	endif
	1310	enddo
	1311
	1312	! for each point within A2, get distance and angle
	1313	! and look where fit with previous table, and set value
	1314
	1315	do ig=1,ngridmx
	1316	val6=dist_pluto(latfi(ig),lonfi(ig),
	1317	& val4pi/180.,val5pi/180.)
	1318	if (val6.lt.val3) then
	1319	! dist = val6
	1320	! get angle:
	1321	angle=atan2(val4-latfi(ig)*180./pi,
	1322	& val5-lonfi(ig)*180./pi)
	1323	! Loop where min
	1324	actualmin=1000.
	1325	do j=1,i
	1326	if ( (array_angle(j).lt.angle+0.52).and.
	1327	& (array_angle(j).gt.angle-0.52).and.
	1328	& (array_dist(j)-val6.lt.actualmin) ) then
	1329	actualmin=array_dist(j)-val6
	1330	igref3=j
	1331	endif
	1332	enddo
	1333	phisfi(ig)=phisfi(array_ind(igref3))
	1334	qsurf(ig,igcm_ch4_ice)=
	1335	& qsurf(array_ind(igref3),igcm_ch4_ice)
	1336	tsurf(ig)=tsurf(array_ind(igref3))
	1337	endif
	1338	enddo
	1339	enddo
	1340	CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
	1341
	1342	c delfrostch4/ delete ch4 frost on a latitudinal band
	1343	c ----------------------------------------------------------
	1344	else if (modif(1:len_trim(modif)) .eq. 'delfrostch4') then
	1345
	1346	write(,) 'removing ch4 latitudinally'
	1347	write(,) 'Choice band : lat1 and lat2 ?'
	1348	read(,,iostat=ierr) val,val2
	1349	write(,) 'Choice band : lon1 and lon2 ?'
	1350	522 read(,,iostat=ierr) val4,val5
	1351	if(ierr.ne.0) goto 522
	1352	write(,) 'Choice amount max below whcih ch4 is removed'
	1353	521 read(,,iostat=ierr) val3
	1354	if(ierr.ne.0) goto 521
	1355
	1356	DO ig=1,ngridmx
	1357	IF ( (latfi(ig)*180./pi.ge.val) .and.
	1358	& (latfi(ig)*180./pi.le.val2) .and.
	1359	& (lonfi(ig)*180./pi.ge.val4) .and.
	1360	& (lonfi(ig)*180./pi.lt.val5) .and.
	1361	& (qsurf(ig,igcm_ch4_ice).lt.val3) ) then
	1362	qsurf(ig,igcm_ch4_ice)=0.
	1363	ENDIF
	1364	ENDDO
	1365
	1366	c modch4 : adding/remove ch4 frost on a latitudinal band
	1367	c ----------------------------------------------------------
	1368	else if (modif(1:len_trim(modif)) .eq. 'modch4') then
	1369
	1370	write(,) 'Adding or removing ch4 latitudinally'
	1371	write(,) 'Choice band : lat1 and lat2 ?'
	1372	read(,,iostat=ierr) val,val2
	1373	write(,) 'Choice band : lon1 and lon2 ?'
	1374	422 read(,,iostat=ierr) val4,val5
	1375	if(ierr.ne.0) goto 422
	1376	write(,) 'Choice multiplicative factor'
	1377	421 read(,,iostat=ierr) val3
	1378	if(ierr.ne.0) goto 421
	1379	write(,) 'Choice additional ch4'
	1380	423 read(,,iostat=ierr) val6
	1381	if(ierr.ne.0) goto 423
	1382
	1383	DO ig=1,ngridmx
	1384	IF ( (latfi(ig)*180./pi.ge.val) .and.
	1385	& (latfi(ig)*180./pi.le.val2) .and.
	1386	& (lonfi(ig)*180./pi.ge.val4) .and.
	1387	& (lonfi(ig)*180./pi.lt.val5)) then
	1388	! & (qsurf(ig,igcm_n2).gt.50)) then
	1389	! & (qsurf(ig,igcm_ch4_ice).lt.10) ) then
	1390	qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)*val3
	1391	qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)+val6
	1392	ENDIF
	1393	ENDDO
	1394
	1395	c modch4n2 : adding/remove ch4 frost on a latitudinal band
	1396	c ----------------------------------------------------------
	1397	else if (modif(1:len_trim(modif)) .eq. 'modch4n2') then
	1398
	1399	write(,) 'Adding or removing ch4 latitudinally'
	1400	write(,) 'Choice band : lat1 and lat2 ?'
	1401	read(,,iostat=ierr) val,val2
	1402	write(,) 'Choice band : lon1 and lon2 ?'
	1403	read(,,iostat=ierr) val4,val5
	1404	write(,) 'Choice range n2 for modif'
	1405	read(,,iostat=ierr) val7,val8
	1406	write(,) 'Choice multiplicative factor ch4'
	1407	read(,,iostat=ierr) val3
	1408	write(,) 'Choice additional ch4'
	1409	read(,,iostat=ierr) val6
	1410
	1411	DO ig=1,ngridmx
	1412	IF ( (latfi(ig)*180./pi.ge.val) .and.
	1413	& (latfi(ig)*180./pi.le.val2) .and.
	1414	& (lonfi(ig)*180./pi.ge.val4) .and.
	1415	& (lonfi(ig)*180./pi.lt.val5) .and.
	1416	& (qsurf(ig,igcm_n2).gt.val7) .and.
	1417	& (qsurf(ig,igcm_n2).lt.val8) ) then
	1418	qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)*val3
	1419	qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)+val6
	1420	ENDIF
	1421	ENDDO
	1422
	1423	c non2noco : if no n2 no co
	1424	c ----------------------------------------------------------
	1425	else if (modif(1:len_trim(modif)) .eq. 'non2noco') then
	1426	DO ig=1,ngridmx
	1427	IF ( (qsurf(ig,igcm_n2).lt.0.001)) then
	1428	qsurf(ig,igcm_co_ice)=0.
	1429	ENDIF
	1430	ENDDO
	1431
	1432	c modco : adding/remove co frost on a latitudinal band
	1433	c ----------------------------------------------------------
	1434	else if (modif(1:len_trim(modif)) .eq. 'modco') then
	1435
	1436	write(,) 'Adding or removing co where N2 is present '
	1437	write(,) 'Choice limit mini n2 pour mettre co?'
	1438	460 read(,,iostat=ierr) val7
	1439	if(ierr.ne.0) goto 460
	1440	write(,) 'Choice band : lat1 and lat2 ?'
	1441	461 read(,,iostat=ierr) val,val2
	1442	if(ierr.ne.0) goto 461
	1443	write(,) 'Choice band : lon1 and lon2 ?'
	1444	462 read(,,iostat=ierr) val4,val5
	1445	if(ierr.ne.0) goto 462
	1446	write(,) 'Choice multiplicative factor'
	1447	463 read(,,iostat=ierr) val3
	1448	if(ierr.ne.0) goto 463
	1449	write(,) 'Choice additional co'
	1450	464 read(,,iostat=ierr) val6
	1451	if(ierr.ne.0) goto 464
	1452
	1453	DO ig=1,ngridmx
	1454	IF ( (latfi(ig)*180./pi.ge.val) .and.
	1455	& (latfi(ig)*180./pi.le.val2) .and.
	1456	& (lonfi(ig)*180./pi.ge.val4) .and.
	1457	& (lonfi(ig)*180./pi.lt.val5) .and.
[3438]	1458	& (qsurf(ig,igcm_n2).gt.val7)) then
[3371]	1459	qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)*val3
	1460	qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)+val6
	1461	ENDIF
	1462	ENDDO
	1463
	1464	c modn2 : adding/remove n2 ice
	1465	c ----------------------------------------------------------
	1466	else if (modif(1:len_trim(modif)) .eq. 'modn2') then
	1467
	1468	write(,) 'Adding or removing n2 '
	1469	write(,) 'Choice band : lat1 and lat2 ?'
	1470	425 read(,,iostat=ierr) val,val2
	1471	if(ierr.ne.0) goto 425
	1472	write(,) 'Choice band : lon1 and lon2 ?'
	1473	426 read(,,iostat=ierr) val4,val5
	1474	if(ierr.ne.0) goto 426
	1475	write(,) 'Choice multiplicative factor'
	1476	427 read(,,iostat=ierr) val3
	1477	if(ierr.ne.0) goto 427
	1478	write(,) 'Choice additional n2'
[3382]	1479	428 read(,,iostat=ierr) val6
[3371]	1480	if(ierr.ne.0) goto 428
	1481
	1482	DO ig=1,ngridmx
	1483	IF ( (latfi(ig)*180./pi.ge.val) .and.
	1484	& (latfi(ig)*180./pi.le.val2) .and.
	1485	& (lonfi(ig)*180./pi.ge.val4) .and.
	1486	& (lonfi(ig)*180./pi.lt.val5) ) then
	1487	c & (qsurf(ig,igcm_n2).lt.50)) then
	1488	qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)*val3
	1489	qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)+val6
	1490	ENDIF
	1491	! IF ((phisfi(ig).gt.-1000.)) then
	1492	! qsurf(ig,igcm_n2)=0.
	1493	! ENDIF
	1494	ENDDO
	1495
	1496	c modcoch4 : adding/remove co and ch4 frost where co without n2
	1497	c ----------------------------------------------------------
	1498	else if (modif(1:len_trim(modif)) .eq. 'modcoch4') then
	1499
	1500	write(,) 'Adding or removing co where N2 is not present '
	1501	write(,) 'Choice band : lat1 and lat2 ?'
	1502	491 read(,,iostat=ierr) val,val2
	1503	if(ierr.ne.0) goto 491
	1504	write(,) 'Choice band : lon1 and lon2 ?'
	1505	492 read(,,iostat=ierr) val4,val5
	1506	if(ierr.ne.0) goto 492
	1507	write(,) 'Choice multiplicative factor'
	1508	493 read(,,iostat=ierr) val3
	1509	if(ierr.ne.0) goto 493
	1510	write(,) 'Choice additional co ch4'
	1511	494 read(,,iostat=ierr) val6
	1512	if(ierr.ne.0) goto 494
	1513
	1514	DO ig=1,ngridmx
	1515	IF ( (latfi(ig)*180./pi.ge.val) .and.
	1516	& (latfi(ig)*180./pi.le.val2) .and.
	1517	& (lonfi(ig)*180./pi.ge.val4) .and.
	1518	& (lonfi(ig)*180./pi.le.val5) .and.
	1519	& (qsurf(ig,igcm_n2).lt.10.)) then
	1520	qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)*val3
	1521	qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)+val6
	1522	qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)*val3
	1523	qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)+val6
	1524	ENDIF
	1525	ENDDO
	1526
	1527	c modn2HD : adding/remove n2 ice for HD runs
	1528	c ----------------------------------------------------------
	1529	else if (modif(1:len_trim(modif)) .eq. 'modn2HD') then
	1530
	1531	write(,) 'Adding or removing n2 '
	1532	write(,) 'Choice band : lat1 and lat2 ?'
	1533	480 read(,,iostat=ierr) val,val1
	1534	if(ierr.ne.0) goto 480
	1535	write(,) 'Choice band : lon1 and lon2 ?'
	1536	481 read(,,iostat=ierr) val2,val3
	1537	if(ierr.ne.0) goto 481
	1538	write(,) 'Choice amount N2 min max where to modify'
	1539	482 read(,,iostat=ierr) val4,val5
	1540	if(ierr.ne.0) goto 482
	1541	write(,) 'Choice phisfi min max where change n2'
	1542	483 read(,,iostat=ierr) val6,val11
	1543	if(ierr.ne.0) goto 483
	1544	write(,) 'Choice multiplicative factor'
	1545	484 read(,,iostat=ierr) val7
	1546	if(ierr.ne.0) goto 484
	1547	write(,) 'Choice additional n2'
	1548	485 read(,,iostat=ierr) val8
	1549	if(ierr.ne.0) goto 485
	1550	! write(,) 'Choice ind lon equivalent for change tsurf tsoil'
	1551	! 486 read(,,iostat=ierr) val9
	1552	! if(ierr.ne.0) goto 486
	1553
	1554	DO ig=1,ngridmx
	1555	IF ( (latfi(ig)*180./pi.ge.val) .and.
	1556	& (latfi(ig)*180./pi.le.val1) .and.
	1557	& (lonfi(ig)*180./pi.ge.val2) .and.
	1558	& (lonfi(ig)*180./pi.lt.val3) .and.
	1559	& (qsurf(ig,igcm_n2).ge.val4) .and.
	1560	& (qsurf(ig,igcm_n2).lt.val5) .and.
	1561	& (phisfi(ig).ge.val6) .and.
	1562	& (phisfi(ig).le.val11) ) then
	1563
	1564	qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)*val7
	1565	qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)+val8
	1566	!qsurf(ig,igcm_ch4_ice)=0.
	1567	qsurf(ig,igcm_co_ice)=0.
	1568	! index of ig
	1569	!if (val9.eq.0.) then
	1570	! iref=2546
	1571	!else
	1572	! val10=modulo((ig-1),96)
	1573	! choice=ig+int(96/2)-val10
	1574	! !choice=int(1+96*(val9-1)+val10)
	1575	!endif
	1576	!tsurf(ig)=tsurf(iref)
	1577	!DO l=1,nsoilmx
	1578	! tsoil(ig,l)=tsoil(iref,l)
	1579	!ENDDO
	1580	tsurf(ig)=tsurf(ig)+4
	1581	DO l=1,nsoilmx
	1582	tsoil(ig,l)=tsoil(ig,l)+4
	1583	ENDDO
	1584	ENDIF
	1585	! IF ((phisfi(ig).gt.-1000.)) then
	1586	! qsurf(ig,igcm_n2)=0.
	1587	! ENDIF
	1588	ENDDO
	1589
	1590	c modn2HD_SP : adding/remove n2 ice for HD runs
	1591	c ----------------------------------------------------------
	1592	else if (modif(1:len_trim(modif)) .eq. 'modn2HD_SP') then
	1593
	1594	write(,) 'Adding or removing n2 '
	1595	write(,) 'Choice band : lat1 and lat2 ?'
	1596	501 read(,,iostat=ierr) val,val1
	1597	if(ierr.ne.0) goto 501
	1598	write(,) 'Choice band : lon1 and lon2 ?'
	1599	502 read(,,iostat=ierr) val2,val3
	1600	if(ierr.ne.0) goto 502
	1601	write(,) 'Choice amount N2 min max where to modify'
	1602	503 read(,,iostat=ierr) val4,val5
	1603	if(ierr.ne.0) goto 503
	1604	write(,) 'Choice phisfi min max where change n2'
[3382]	1605	504 read(,,iostat=ierr) val6,val11
[3371]	1606	if(ierr.ne.0) goto 504
	1607	write(,) 'Choice multiplicative factor'
	1608	505 read(,,iostat=ierr) val7
	1609	if(ierr.ne.0) goto 505
	1610	write(,) 'Choice additional n2'
[3382]	1611	506 read(,,iostat=ierr) val8
[3371]	1612	if(ierr.ne.0) goto 506
	1613	write(,) 'Choice ind lon equivalent for change tsurf tsoil'
[3382]	1614	507 read(,,iostat=ierr) iref
[3371]	1615	if(ierr.ne.0) goto 507
	1616
	1617	DO ig=1,ngridmx
	1618	IF ( (latfi(ig)*180./pi.ge.val) .and.
	1619	& (latfi(ig)*180./pi.le.val1) .and.
	1620	& (lonfi(ig)*180./pi.ge.val2) .and.
	1621	& (lonfi(ig)*180./pi.lt.val3) .and.
	1622	& (qsurf(ig,igcm_n2).ge.val4) .and.
	1623	& (qsurf(ig,igcm_n2).lt.val5) .and.
	1624	& (phisfi(ig).ge.val6) .and.
	1625	& (phisfi(ig).le.val11) ) then
	1626
	1627	qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)*val7
	1628	qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)+val8
	1629	qsurf(ig,igcm_ch4_ice)=0.
	1630	qsurf(ig,igcm_co_ice)=0.
	1631	! index of ig
	1632	!if (val9.eq.0.) then
	1633	! iref=2546
	1634	!else
	1635	!val10=modulo((ig-1),96)
	1636	!choice=ig+int(96/2)-val10
	1637	!choice=int(1+96*(val9-1)+val10)
	1638	if (iref.ge.1) then
	1639	tsurf(ig)=tsurf(iref)
	1640	DO l=1,nsoilmx
	1641	tsoil(ig,l)=tsoil(iref,l)
	1642	ENDDO
	1643	else if (iref.eq.0) then
	1644	choice=int(ig/96.)*96+84
	1645	print*, 'choice=',choice
	1646	tsurf(ig)=tsurf(int(choice))
	1647	DO l=1,nsoilmx
	1648	tsoil(ig,l)=tsoil(int(choice),l)
	1649	ENDDO
	1650	endif
	1651	ENDIF
	1652	ENDDO
	1653
	1654	c modn2_topo : adding/remove n2 ice
	1655	c ----------------------------------------------------------
	1656	else if (modif(1:len_trim(modif)) .eq. 'modn2_topo') then
	1657
	1658	write(,) 'Adding or removing n2 '
	1659	write(,) 'Choice band : lat1 and lat2 ?'
	1660	441 read(,,iostat=ierr) val,val2
	1661	if(ierr.ne.0) goto 441
	1662	write(,) 'Choice band : lon1 and lon2 ?'
	1663	442 read(,,iostat=ierr) val4,val5
	1664	if(ierr.ne.0) goto 442
	1665	write(,) 'Choice topo 1 et 2 (phi) where change is active'
	1666	443 read(,,iostat=ierr) val7,val8
	1667	if(ierr.ne.0) goto 443
	1668	write(,) 'Choice multiplicative factor'
	1669	444 read(,,iostat=ierr) val3
	1670	if(ierr.ne.0) goto 444
	1671	write(,) 'Choice additional n2'
[3382]	1672	445 read(,,iostat=ierr) val6
[3371]	1673	if(ierr.ne.0) goto 445
	1674
	1675	DO ig=1,ngridmx
	1676	IF ( (latfi(ig)*180./pi.ge.val) .and.
	1677	& (latfi(ig)*180./pi.le.val2) .and.
	1678	& (lonfi(ig)*180./pi.ge.val4) .and.
	1679	& (lonfi(ig)*180./pi.lt.val5) .and.
	1680	& (phisfi(ig).le.val8) .and.
	1681	& (phisfi(ig).ge.val7) ) then
	1682	qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)*val3
	1683	qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)+val6
	1684	ENDIF
	1685	ENDDO
	1686
	1687	c modwheren2 : adding/remove n2 ice where already n2
	1688	c ----------------------------------------------------------
	1689	else if (modif(1:len_trim(modif)) .eq. 'modwheren2') then
	1690
	1691	write(,) 'Choice multiplicative factor'
	1692	471 read(,,iostat=ierr) val3
	1693	if(ierr.ne.0) goto 471
	1694	write(,) 'Choice additional n2'
[3382]	1695	472 read(,,iostat=ierr) val6
[3371]	1696	if(ierr.ne.0) goto 472
	1697
	1698	DO ig=1,ngridmx
	1699	IF ( qsurf(ig,igcm_n2).gt.0.001 ) then
	1700	qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)*val3
	1701	qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)+val6
	1702	ENDIF
	1703	ENDDO
	1704
	1705	c checktsoil : changing tsoil if no n2
	1706	c ----------------------------------------------------------
	1707	else if (modif(1:len_trim(modif)) .eq. 'checktsoil') then
	1708
	1709	write(,) 'selecting area where tsoil to be check '
	1710	write(,) 'Choice band : lat1 and lat2 ?'
	1711	511 read(,,iostat=ierr) val,val1
	1712	if(ierr.ne.0) goto 511
	1713	write(,) 'Choice band : lon1 and lon2 ?'
	1714	512 read(,,iostat=ierr) val2,val3
	1715	if(ierr.ne.0) goto 512
	1716	write(,) 'Choice temp min max in which change is made'
	1717	513 read(,,iostat=ierr) val4,val5
	1718	if(ierr.ne.0) goto 513
	1719	write(,) 'Choice phisfi min max where change n2'
[3382]	1720	514 read(,,iostat=ierr) val6,val11
[3371]	1721	if(ierr.ne.0) goto 514
	1722
	1723	DO ig=1,ngridmx
	1724	IF ( (latfi(ig)*180./pi.ge.val) .and.
	1725	& (latfi(ig)*180./pi.le.val1) .and.
	1726	& (lonfi(ig)*180./pi.ge.val2) .and.
	1727	& (lonfi(ig)*180./pi.le.val3) .and.
	1728	& (((tsurf(ig).ge.val4) .and.
	1729	& (tsurf(ig).le.val5)) .or.
	1730	& ((tsoil(ig,nsoilmx).ge.val4) .and.
	1731	& (tsoil(ig,nsoilmx).le.val5))) .and.
	1732	& (phisfi(ig).ge.val6) .and.
	1733	& (phisfi(ig).le.val11) .and.
[3390]	1734	& (qsurf(ig,igcm_n2).gt.0.001) ) then
[3371]	1735
	1736	! DO i=1,ngridmx
	1737	! IF ( (latfi(i).eq.latfi(ig)) .and.
	1738	! & (lonfi(i).eq.0.) ) then
	1739	!
[3390]	1740	tsurf(ig)=34.7
[3371]	1741	!qsurf(ig,igcm_ch4_ice)=qsurf(i,igcm_ch4_ice)
	1742	!
	1743	DO l=1,nsoilmx
[3390]	1744	tsoil(ig,l)=34.7 !tsoil(i,l)
[3371]	1745	ENDDO
	1746	!ENDIF
	1747	!ENDDO
	1748	ENDIF
	1749	ENDDO
	1750
	1751	c asymtriton : changing ice, tsurf and tsoil on triton
	1752	c ----------------------------------------------------------
	1753	else if (modif(1:len_trim(modif)) .eq. 'asymtriton') then
	1754
	1755	write(,) 'selecting area where tsoil to be check '
	1756	write(,) 'Choice center latitude of sinuisoid distrib?'
	1757	531 read(,,iostat=ierr) val1
	1758	if(ierr.ne.0) goto 531
	1759	write(,) 'Choice maginutde in latitude?'
	1760	532 read(,,iostat=ierr) val2
	1761	if(ierr.ne.0) goto 532
	1762	write(,) 'Choice center longitude?'
	1763	533 read(,,iostat=ierr) val3
	1764	if(ierr.ne.0) goto 533
	1765	! write(,) 'iip1,jjp1=',iip1,jjp1,ngridmx
	1766
	1767	DO ig=1,ngridmx
	1768	! Latitude of the sinusoid:
	1769	val11=val1+val2cos(lonfi(ig)1.570796332/pi-val3pi/180.)
	1770	! If above line and ice: remove ice
	1771	IF ( (latfi(ig)*180./pi.ge.val11) .and.
	1772	& (latfi(ig)*180./pi.le.val1+val2) .and.
	1773	& (qsurf(ig,igcm_n2).gt.0.) ) then
	1774	! Look for same longitude point where no ice
	1775	val5=1.
	1776	jref=ig
	1777	! 1
	1778	! ... iip1 ... x (jjp1-2) 32 x 23
	1779	! 1
	1780	do while (val5.ge.1..and.jref.gt.iip1)
	1781	! northward point
	1782	jref=jref-iip1+1
	1783	! ice at the point
	1784	val5=qsurf(jref,igcm_n2)
	1785	! write(,) jref,
	1786	! & latfi(jref)180./pi,lonfi(jref)180/pi,val5
	1787	enddo
	1788	if (val5.ge.1) write(,) 'NO POINT FOUND WITH NO ICE'
	1789
	1790	! Copy point in the selected area
	1791	tsurf(ig)=tsurf(jref)
	1792	qsurf(ig,igcm_n2)=qsurf(jref,igcm_n2)
	1793	qsurf(ig,igcm_ch4_ice)=qsurf(jref,igcm_ch4_ice)
	1794	qsurf(ig,igcm_co_ice)=qsurf(jref,igcm_co_ice)
	1795	DO l=1,nsoilmx
	1796	tsoil(ig,l)=tsoil(jref,l)
	1797	ENDDO
	1798	ENDIF
	1799	! If below line and no ice: add ice
	1800	IF ( (latfi(ig)*180./pi.le.val11) .and.
	1801	& (latfi(ig)*180./pi.le.val1+val2) .and.
	1802	& (qsurf(ig,igcm_n2).eq.0.) ) then
	1803	! Look for same longitude point where ice
	1804	val5=1.
	1805	jref=ig
	1806	do while (val5.le.1.and.jref.lt.ngridmx-iip1)
	1807	! southward point
	1808	jref=jref+iip1-1
	1809	! ice at the point
	1810	val5=qsurf(jref,igcm_n2)
	1811	write(,) jref,
	1812	& latfi(jref)180./pi,lonfi(jref)180/pi,val5
	1813	enddo
	1814	if (val5.le.1) write(,) 'NO POINT FOUND WITH ICE'
	1815
	1816	! Copy point in the selected area
	1817	tsurf(ig)=tsurf(jref)
	1818	qsurf(ig,igcm_n2)=qsurf(jref,igcm_n2)
	1819	qsurf(ig,igcm_ch4_ice)=qsurf(jref,igcm_ch4_ice)
	1820	qsurf(ig,igcm_co_ice)=qsurf(jref,igcm_co_ice)
	1821	DO l=1,nsoilmx
	1822	tsoil(ig,l)=tsoil(jref,l)
	1823	ENDDO
	1824	ENDIF
	1825
	1826	ENDDO
	1827
	1828	c source_ch4 : adding source ch4
	1829	c ----------------------------------------------------------
	1830	else if (modif(1:len_trim(modif)) .eq. 'source_ch4') then
	1831
	1832	write(,) 'Adding ch4 ice latitudinally '
	1833	write(,) 'Choice : lat1 and lat2 ?'
	1834	433 read(,,iostat=ierr) val,val2
	1835	if(ierr.ne.0) goto 433
	1836	write(,) 'Choice : lon1 and lon2 ?'
	1837	434 read(,,iostat=ierr) val3,val4
	1838	if(ierr.ne.0) goto 434
	1839	write(,) 'Choice amount (kg/m2)'
	1840	435 read(,,iostat=ierr) val5
	1841	if(ierr.ne.0) goto 435
	1842
	1843	DO ig=1,ngridmx
	1844	IF ( (latfi(ig)*180./pi.ge.val) .and.
	1845	& (latfi(ig)*180./pi.le.val2) .and.
	1846	& (lonfi(ig)*180./pi.ge.val3) .and.
	1847	& (lonfi(ig)*180./pi.lt.val4) ) then
	1848	qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)+val5
	1849	ENDIF
	1850	ENDDO
	1851
	1852	c source_co : adding source co
	1853	c ----------------------------------------------------------
	1854	else if (modif(1:len_trim(modif)) .eq. 'source_co') then
	1855
	1856	write(,) 'Adding co ice latitudinally '
	1857	write(,) 'Choice : lat1 and lat2 ?'
	1858	436 read(,,iostat=ierr) val,val2
	1859	if(ierr.ne.0) goto 436
	1860	write(,) 'Choice : lon1 and lon2 ?'
	1861	437 read(,,iostat=ierr) val3,val4
	1862	if(ierr.ne.0) goto 437
	1863	write(,) 'Choice amount (kg/m2)'
	1864	438 read(,,iostat=ierr) val5
	1865	if(ierr.ne.0) goto 438
	1866
	1867	DO ig=1,ngridmx
	1868	IF ( (latfi(ig)*180./pi.ge.val) .and.
	1869	& (latfi(ig)*180./pi.le.val2) .and.
	1870	& (lonfi(ig)*180./pi.ge.val3) .and.
	1871	& (lonfi(ig)*180./pi.lt.val4) ) then
	1872	qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)+val5
	1873	ENDIF
	1874	ENDDO
	1875
	1876	! therm_ini_s: (re)-set soil thermal inertia to reference surface values
	1877	! ----------------------------------------------------------------------
	1878	else if (modif(1:len_trim(modif)).eq.'therm_ini_s') then
	1879	! write(,)"surfithfi(1):",surfithfi(1)
	1880	do isoil=1,nsoilmx
	1881	inertiedat(1:ngridmx,isoil)=surfithfi(1:ngridmx)
	1882	enddo
	1883	write(,)'OK: Soil thermal inertia has been reset to referenc
[3382]	1884	& e surface values'
[3371]	1885	write(,)"inertiedat(1,1):",inertiedat(1,1)
	1886	ithfi(:,:)=inertiedat(:,:)
	1887	! recast ithfi() onto ith()
	1888	call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
	1889
	1890	! inert3d: set soil thermal inertia to specific uniform value
	1891	! ----------------------------------------------------------------------
	1892	else if (modif(1:len_trim(modif)).eq.'inert3d') then
	1893	write(,) 'Actual value of surf Thermal inertia at ig=1: ',
	1894	& inertiedat(1,1), ' SI'
	1895	write(,) 'Value of Thermal inertia (SI) ? '
	1896	read(,) val
	1897	do isoil=1,nsoilmx
	1898	do ig=1,ngridmx
	1899	inertiedat(ig,isoil)=val
	1900	enddo
	1901	enddo
	1902	write(,)'OK: Soil TI set to ',val,' SI everywhere'
	1903	ithfi(:,:)=inertiedat(:,:)
	1904	! recast ithfi() onto ith()
	1905	call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
	1906
	1907	! subsoilice_n: Put deep ice layer in northern hemisphere soil
	1908	! ------------------------------------------------------------
[3382]	1909	else if (modif(1:len_trim(modif)).eq.'subsoilice_n') then
[3371]	1910
	1911	write(,)'From which latitude (in deg.), up to the north pole,
[3382]	1912	& should we put subterranean ice?'
[3228]	1913	ierr=1
	1914	do while (ierr.ne.0)
[3371]	1915	read(,,iostat=ierr) val
	1916	if (ierr.eq.0) then ! got a value
	1917	! do a sanity check
	1918	if((val.lt.0.).or.(val.gt.90)) then
	1919	write(,)'Latitude should be between 0 and 90 deg. !!!'
	1920	ierr=1
	1921	else ! find corresponding jref (nearest latitude)
	1922	! note: rlatu(:) contains decreasing values of latitude
	1923	! starting from PI/2 to -PI/2
	1924	do j=1,jjp1
	1925	if ((rlatu(j)*180./pi.ge.val).and.
	1926	& (rlatu(j+1)*180./pi.le.val)) then
	1927	! find which grid point is nearest to val:
	1928	if (abs(rlatu(j)*180./pi-val).le.
	1929	& abs((rlatu(j+1)*180./pi-val))) then
	1930	jref=j
	1931	else
	1932	jref=j+1
	1933	endif
	1934
	1935	write(,)'Will use nearest grid latitude which is:',
	1936	& rlatu(jref)*180./pi
	1937	endif
	1938	enddo ! of do j=1,jjp1
	1939	endif ! of if((val.lt.0.).or.(val.gt.90))
	1940	endif !of if (ierr.eq.0)
	1941	enddo ! of do while
	1942
	1943	! Build layers() (as in soil_settings.F)
	1944	val2=sqrt(mlayer(0)*mlayer(1))
	1945	val3=mlayer(1)/mlayer(0)
	1946	do isoil=1,nsoilmx
	1947	layer(isoil)=val2(val3*(isoil-1))
	1948	enddo
	1949
	1950	write(,)'At which depth (in m.) does the ice layer begin?'
	1951	write(,)'(currently, the deepest soil layer extends down to:'
	1952	& ,layer(nsoilmx),')'
	1953	ierr=1
	1954	do while (ierr.ne.0)
[3228]	1955	read(,,iostat=ierr) val2
[3371]	1956	! write(,)'val2:',val2,'ierr=',ierr
[3228]	1957	if (ierr.eq.0) then ! got a value, but do a sanity check
	1958	if(val2.gt.layer(nsoilmx)) then
	1959	write(,)'Depth should be less than ',layer(nsoilmx)
	1960	ierr=1
	1961	endif
	1962	if(val2.lt.layer(1)) then
	1963	write(,)'Depth should be more than ',layer(1)
	1964	ierr=1
	1965	endif
	1966	endif
	1967	enddo ! of do while
[3371]	1968
	1969	! find the reference index iref the depth corresponds to
	1970	! if (val2.lt.layer(1)) then
	1971	! iref=1
	1972	! else
[3228]	1973	do isoil=1,nsoilmx-1
[3371]	1974	if((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1)))
[3228]	1975	& then
[3371]	1976	iref=isoil
[3228]	1977	exit
[3371]	1978	endif
[3228]	1979	enddo
	1980
[3371]	1981	! thermal inertia of the ice:
	1982	ierr=1
	1983	do while (ierr.ne.0)
	1984	write(,)'What is the value of subterranean ice thermal inert
[3382]	1985	& ia? (e.g.: 2000)'
[3371]	1986	read(,,iostat=ierr)iceith
	1987	enddo ! of do while
[3228]	1988
[3371]	1989	! recast ithfi() onto ith()
	1990	call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
	1991
	1992	do j=1,jref
	1993	! write(,)'j:',j,'rlatu(j)180./pi:',rlatu(j)180./pi
	1994	do i=1,iip1 ! loop on longitudes
	1995	! Build "equivalent" thermal inertia for the mixed layer
	1996	ith(i,j,iref+1)=sqrt((layer(iref+1)-layer(iref))/
	1997	& (((val2-layer(iref))/(ith(i,j,iref)**2))+
	1998	& ((layer(iref+1)-val2)/(iceith)**2)))
	1999	! Set thermal inertia of lower layers
	2000	do isoil=iref+2,nsoilmx
	2001	ith(i,j,isoil)=iceith ! ice
	2002	enddo
	2003	enddo ! of do i=1,iip1
	2004	enddo ! of do j=1,jjp1
	2005
	2006	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
[3228]	2007
[3371]	2008	! subsoilice_s: Put deep ice layer in southern hemisphere soil
	2009	! ------------------------------------------------------------
	2010	else if (modif(1:len_trim(modif)).eq.'subsoilice_s') then
[3228]	2011
[3371]	2012	write(,)'From which latitude (in deg.), down to the south pol
[3382]	2013	& e, should we put subterranean ice?'
[3371]	2014	ierr=1
	2015	do while (ierr.ne.0)
	2016	read(,,iostat=ierr) val
	2017	if (ierr.eq.0) then ! got a value
	2018	! do a sanity check
	2019	if((val.gt.0.).or.(val.lt.-90)) then
	2020	write(,)'Latitude should be between 0 and -90 deg. !!!'
	2021	ierr=1
	2022	else ! find corresponding jref (nearest latitude)
	2023	! note: rlatu(:) contains decreasing values of latitude
	2024	! starting from PI/2 to -PI/2
	2025	do j=1,jjp1
	2026	if ((rlatu(j)*180./pi.ge.val).and.
	2027	& (rlatu(j+1)*180./pi.le.val)) then
	2028	! find which grid point is nearest to val:
	2029	if (abs(rlatu(j)*180./pi-val).le.
	2030	& abs((rlatu(j+1)*180./pi-val))) then
	2031	jref=j
	2032	else
	2033	jref=j+1
	2034	endif
	2035
	2036	write(,)'Will use nearest grid latitude which is:',
	2037	& rlatu(jref)*180./pi
	2038	endif
	2039	enddo ! of do j=1,jjp1
	2040	endif ! of if((val.lt.0.).or.(val.gt.90))
	2041	endif !of if (ierr.eq.0)
	2042	enddo ! of do while
[3228]	2043
[3371]	2044	! Build layers() (as in soil_settings.F)
	2045	val2=sqrt(mlayer(0)*mlayer(1))
	2046	val3=mlayer(1)/mlayer(0)
	2047	do isoil=1,nsoilmx
	2048	layer(isoil)=val2(val3*(isoil-1))
	2049	enddo
[3228]	2050
[3371]	2051	write(,)'At which depth (in m.) does the ice layer begin?'
	2052	write(,)'(currently, the deepest soil layer extends down to:'
	2053	& ,layer(nsoilmx),')'
	2054	ierr=1
	2055	do while (ierr.ne.0)
	2056	read(,,iostat=ierr) val2
	2057	! write(,)'val2:',val2,'ierr=',ierr
	2058	if (ierr.eq.0) then ! got a value, but do a sanity check
	2059	if(val2.gt.layer(nsoilmx)) then
	2060	write(,)'Depth should be less than ',layer(nsoilmx)
	2061	ierr=1
	2062	endif
	2063	if(val2.lt.layer(1)) then
	2064	write(,)'Depth should be more than ',layer(1)
	2065	ierr=1
	2066	endif
	2067	endif
	2068	enddo ! of do while
	2069
	2070	! find the reference index iref the depth corresponds to
	2071	do isoil=1,nsoilmx-1
	2072	if((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1)))
	2073	& then
	2074	iref=isoil
	2075	exit
	2076	endif
	2077	enddo
[3382]	2078
[3371]	2079	! thermal inertia of the ice:
	2080	ierr=1
	2081	do while (ierr.ne.0)
	2082	write(,)'What is the value of subterranean ice thermal inert
[3382]	2083	& ia? (e.g.: 2000)'
[3371]	2084	read(,,iostat=ierr)iceith
	2085	enddo ! of do while
	2086
	2087	! recast ithfi() onto ith()
	2088	call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
	2089
	2090	do j=jref,jjp1
	2091	! write(,)'j:',j,'rlatu(j)180./pi:',rlatu(j)180./pi
	2092	do i=1,iip1 ! loop on longitudes
	2093	! Build "equivalent" thermal inertia for the mixed layer
	2094	ith(i,j,iref+1)=sqrt((layer(iref+1)-layer(iref))/
	2095	& (((val2-layer(iref))/(ith(i,j,iref)**2))+
	2096	& ((layer(iref+1)-val2)/(iceith)**2)))
	2097	! Set thermal inertia of lower layers
	2098	do isoil=iref+2,nsoilmx
	2099	ith(i,j,isoil)=iceith ! ice
	2100	enddo
	2101	enddo ! of do i=1,iip1
	2102	enddo ! of do j=jref,jjp1
[3228]	2103
[3371]	2104	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
[3228]	2105
[3371]	2106	c ----------------------------------------------------------
	2107	c 'reservoir_SP' : increase or decrase ices reservoir in SP by factor
	2108	c ----------------------------------------------------------
	2109	else if (modif(1:len_trim(modif)).eq.'reservoir_SP') then
	2110	write(,) "Increase/Decrease N2 reservoir by factor :"
	2111	read(,) val7
	2112	write(,) "Increase/Decrease CH4 reservoir by factor :"
	2113	read(,) val8
	2114	write(,) "Increase/Decrease CO reservoir by factor :"
	2115	read(,) val9
[3228]	2116
[3371]	2117	! Definition SP:
[3382]	2118	val3=-33. !lat1
	2119	val4=50. !lat2
	2120	val5=140. ! lon1
	2121	val6=-155. ! lon2
	2122	choice=-50. ! min phisfi
	2123	write(,) 'def SP :'
	2124	write(,) 'lat :',val3,val4
	2125	write(,) 'lon :',val5,'180 / -180',val6
	2126	write(,) 'choice phisfi min :',choice
	2127
	2128	DO ig=1,ngridmx
	2129
	2130	IF ( (latfi(ig)*180./pi.ge.val3) .and.
	2131	& (latfi(ig)*180./pi.le.val4) .and.
	2132	& ( ((lonfi(ig)*180./pi.ge.-180.) .and.
	2133	& (lonfi(ig)*180./pi.le.val6)) .or.
	2134	& ((lonfi(ig)*180./pi.ge.val5) .and.
	2135	& (lonfi(ig)*180./pi.le.180.))) ) then
	2136
	2137	IF ((phisfi(ig).le.choice)) then !.and.
	2138	c & (qsurf(ig,igcm_n2).ge.50)) then
	2139
	2140	qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)*val7
	2141	qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)*val8
	2142	qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)*val9
	2143	ENDIF
	2144	ENDIF
	2145	ENDDO
	2146
	2147	c subsoil_SP : choice of thermal inertia values for SP
	2148	c ----------------------------------------------------------------
	2149	else if (modif(1:len_trim(modif)) .eq. 'subsoil_SP') then
	2150
	2151	write(,) 'New value for subsoil thermal inertia in SP ?'
	2152	510 read(,,iostat=ierr) ith_bb
	2153	if(ierr.ne.0) goto 510
	2154	write(,) 'thermal inertia (new value):',ith_bb
	2155
	2156	write(,)'At which depth (in m.) does the ice layer begin?'
	2157	write(,)'(currently, the deepest soil layer extends down to:'
	2158	& ,layer(1),' - ',layer(nsoilmx),')'
	2159	write(,)'write 0 for uniform value for all subsurf levels?'
	2160	ierr=1
	2161	do while (ierr.ne.0)
	2162	read(,,iostat=ierr) val2
	2163	write(,)'val2:',val2,'ierr=',ierr
	2164	if (ierr.eq.0) then ! got a value, but do a sanity check
	2165	if(val2.gt.layer(nsoilmx)) then
	2166	write(,)'Depth should be less than ',layer(nsoilmx)
	2167	ierr=1
	2168	endif
	2169	if(val2.lt.layer(1)) then
	2170	if(val2.eq.0) then
	2171	write(,)'Thermal inertia set for all subsurface layers'
	2172	ierr=0
	2173	else
	2174	write(,)'Depth should be more than ',layer(1)
	2175	ierr=1
	2176	endif
	2177	endif
	2178	endif
	2179	enddo ! of do while
	2180
	2181	! find the reference index iref the depth corresponds to
	2182	if(val2.eq.0) then
	2183	iref=1
	2184	write(,)'Level selected is first level: ',layer(iref),' m'
	2185	else
	2186	do isoil=1,nsoilmx-1
	2187	if((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1)))
	2188	& then
	2189	iref=isoil
	2190	write(,)'Level selected : ',layer(isoil),' m'
	2191	exit
	2192	endif
	2193	enddo
	2194	endif
	2195
	2196	! Definition SP:
	2197	val3=-50. !lat1
	2198	val4=60. !lat2
	2199	val5=130. ! lon1
	2200	val6=-140. ! lon2
	2201	choice=-100. ! min phisfi
[3371]	2202	write(,) 'def SP :'
	2203	write(,) 'lat :',val3,val4
	2204	write(,) 'lon :',val5,'180 / -180',val6
	2205	write(,) 'choice phisfi min :',choice
[3228]	2206
[3371]	2207	DO ig=1,ngridmx
[3228]	2208
[3371]	2209	IF ( (latfi(ig)*180./pi.ge.val3) .and.
	2210	& (latfi(ig)*180./pi.le.val4) .and.
	2211	& ( ((lonfi(ig)*180./pi.ge.-180.) .and.
	2212	& (lonfi(ig)*180./pi.le.val6)) .or.
	2213	& ((lonfi(ig)*180./pi.ge.val5) .and.
	2214	& (lonfi(ig)*180./pi.le.180.))) ) then
[3228]	2215
[3382]	2216	IF ((phisfi(ig).le.choice) .and.
	2217	& (qsurf(ig,igcm_n2).ge.1000)) then
	2218	DO j=iref,nsoilmx
	2219	ithfi(ig,j)=ith_bb
	2220	ENDDO
	2221	ENDIF
	2222	ENDIF
	2223	ENDDO
[3228]	2224
[3382]	2225	CALL gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
	2226
	2227	c topo_sp : change topo of SP
	2228	c -----------------------------------------------------
	2229	else if (modif(1:len_trim(modif)) .eq. 'topo_sp') then
	2230
	2231	! def SP:
	2232	val2=-33. !lat1
	2233	val3=50. !lat2
	2234	val4=140. ! lon1
	2235	val5=-155. ! lon2
	2236	write(,) 'def SP :'
	2237	write(,) 'lat :',val2,val3
	2238	write(,) 'lon :',val4,'180 / -180',val5
	2239	write(,) 'choice phisfi min (ex: -100):'
	2240	606 read(,,iostat=ierr) val6
	2241	if(ierr.ne.0) goto 606
	2242	write(,) 'choice coefficient for phisfi (ex: 2):'
	2243	605 read(,,iostat=ierr) choice
	2244	if(ierr.ne.0) goto 605
	2245
	2246	write(,) 're scale the pressure :'
	2247	r = 8.314511E+0 *1000.E+0/mugaz
	2248	temp=40.
	2249	write(,) 'r, sale height temperature =',r,temp
	2250
	2251	do j=1,jjp1
	2252	do i=1,iip1
	2253	phisprev= phis(i,j)
	2254	IF ( (rlatu(j)*180./pi.ge.val2) .and.
	2255	& (rlatu(j)*180./pi.le.val3) .and.
	2256	& ( ((rlonv(i)*180./pi.ge.-180.) .and.
	2257	& (rlonv(i)*180./pi.le.val5)) .or.
	2258	& ((rlonv(i)*180./pi.ge.val4) .and.
	2259	& (rlonv(i)*180./pi.le.180.))) ) then
	2260
	2261	IF (phis(i,j).le.val6) then
	2262	phis(i,j)=phis(i,j)*choice
	2263	ps(i,j) = ps(i,j) *
	2264	. exp((phisprev-phis(i,j))/(temp * r))
	2265	ENDIF
	2266	ENDIF
	2267	end do
	2268	end do
	2269	c periodicity of surface ps in longitude
	2270	do j=1,jjp1
	2271	ps(1,j) = ps(iip1,j)
	2272	end do
	2273	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
	2274
	2275	c fill_sp inv: fill inverted SP with N2 ice and adjust phisfi (flat SP)
	2276	c -----------------------------------------------------
	2277	else if (modif(1:len_trim(modif)) .eq. 'fill_sp_inv') then
	2278
	2279	! def SP:
	2280	val2=-33. !lat1
	2281	val3=50. !lat2
	2282	val4=-40. ! lon1
	2283	val5=25. ! lon2
	2284	write(,) 'def SP :'
	2285	write(,) 'lat :',val2,val3
	2286	write(,) 'lon :',val4,val5
	2287	write(,) 'choice phisfi ref SP (ex: -800):'
	2288	706 read(,,iostat=ierr) val6
	2289	if(ierr.ne.0) goto 706
	2290
	2291	write(,) 're scale the pressure :'
	2292	r = 8.314511E+0 *1000.E+0/mugaz
	2293	temp=40.
	2294	write(,) 'r, sale height temperature =',r,temp,g
	2295	qsurf=0.
	2296
	2297	write(,) 'latfi=',latfi
	2298	!CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
	2299	write(,) 'phis=',phis
	2300	write(,) 'phisfi=',phisfi
	2301	!CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
	2302	phisold = phis
	2303	write(,) 'phisold=',phisold
	2304	do ig=1,ngridmx
	2305
	2306	write(,) 'lat=',latfi(ig)*180./pi
	2307	write(,) 'lon=',lonfi(ig)*180./pi
	2308	write(,) 'phisfi=',phisfi(ig)
	2309	IF ( (latfi(ig)*180./pi.ge.val2) .and.
	2310	& (latfi(ig)*180./pi.le.val3) .and.
	2311	& (lonfi(ig)*180./pi.ge.val4) .and.
	2312	& (lonfi(ig)*180./pi.le.val5) ) then
	2313	write(,) 'hello SP',phisfi(ig),ig
	2314	IF (phisfi(ig).le.val6) then
	2315	qsurf(ig,igcm_n2)=(val6-phisfi(ig))/g*1000.
	2316	phisfi(ig)=val6
	2317	write(,) 'changes',phisfi(ig),qsurf(ig,igcm_n2)
[3371]	2318	ENDIF
[3382]	2319	ENDIF
	2320	end do
	2321	c update new phis and ps
	2322	CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
	2323	do j=1,jjp1
	2324	do i=1,iip1
	2325	ps(i,j) = ps(i,j) *
	2326	. exp((phisold(i,j)-phis(i,j))/(temp * r))
	2327	enddo
	2328	enddo
	2329	c periodicity of surface ps in longitude
	2330	do j=1,jjp1
	2331	ps(1,j) = ps(iip1,j)
	2332	end do
	2333
	2334	c adjust phisfi according to the amount of N2 ice
	2335	c -----------------------------------------------------
	2336	else if (modif(1:len_trim(modif)) .eq. 'adjustphi') then
	2337
	2338	phisold = phis
	2339	do ig=1,ngridmx
	2340	phisfi(ig)=phisfi(ig)+qsurf(ig,igcm_n2)*g/1000.
	2341	end do
	2342	c update new phis and ps
	2343	CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
	2344	r = 8.314511E+0 *1000.E+0/mugaz
	2345	temp=37.
	2346	do j=1,jjp1
	2347	do i=1,iip1
	2348	ps(i,j) = ps(i,j) *
	2349	. exp((phisold(i,j)-phis(i,j))/(temp * r))
	2350	enddo
	2351	enddo
	2352	c periodicity of surface ps in longitude
	2353	do j=1,jjp1
	2354	ps(1,j) = ps(iip1,j)
	2355	end do
	2356
	2357	c correct phisfi
	2358	c -----------------------------------------------------
	2359	else if (modif(1:len_trim(modif)) .eq. 'correctphi') then
	2360
	2361	write(,) 'choice latitudes:'
	2362	537 read(,,iostat=ierr) val1,val2
	2363	if(ierr.ne.0) goto 537
	2364	write(,) 'choice longitudes:'
	2365	538 read(,,iostat=ierr) val3,val4
	2366	if(ierr.ne.0) goto 538
	2367	write(,) 'choice phi min max'
	2368	539 read(,,iostat=ierr) val5,val6
	2369	if(ierr.ne.0) goto 539
	2370	write(,) 'choice factor phi'
	2371	535 read(,,iostat=ierr) val8
	2372	if(ierr.ne.0) goto 535
	2373	write(,) 'choice add phi'
	2374	536 read(,,iostat=ierr) val7
	2375	if(ierr.ne.0) goto 536
	2376
	2377	do ig=1,ngridmx
	2378	IF ( (latfi(ig)*180./pi.ge.val1) .and.
	2379	& (latfi(ig)*180./pi.le.val2) .and.
	2380	& (lonfi(ig)*180./pi.ge.val3) .and.
	2381	& (lonfi(ig)*180./pi.le.val4) ) then
	2382
	2383	IF ((phisfi(ig).le.val6).and.
	2384	& (phisfi(ig).ge.val5)) then
	2385
	2386	phisfi(ig)=phisfi(ig)*val8
	2387	phisfi(ig)=phisfi(ig)+val7
	2388	ENDIF
	2389
	2390	ENDIF
	2391	enddo
	2392	phisold = phis
	2393	c update new phis and ps
	2394	CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
	2395	r = 8.314511E+0 *1000.E+0/mugaz
	2396	temp=37.
	2397	do j=1,jjp1
	2398	do i=1,iip1
	2399	ps(i,j) = ps(i,j) *
	2400	. exp((phisold(i,j)-phis(i,j))/(temp * r))
	2401	enddo
	2402	enddo
	2403	c periodicity of surface ps in longitude
	2404	do j=1,jjp1
	2405	do i=1,iip1
	2406	ps(i,j) = ps(i,j) *
	2407	. exp((phisold(i,j)-phis(i,j))/(temp * r))
	2408	enddo
	2409	enddo
	2410	c periodicity of surface ps in longitude
	2411	do j=1,jjp1
	2412	ps(1,j) = ps(iip1,j)
	2413	end do
	2414	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
	2415
	2416	c correct phisfi
	2417	c -----------------------------------------------------
	2418	else if (modif(1:len_trim(modif)) .eq. 'correctps') then
	2419
	2420	phisold = phis
	2421	c update new phis and ps
	2422	CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
	2423	r = 8.314511E+0 *1000.E+0/mugaz
	2424	temp=37.
	2425	do j=1,jjp1
	2426	do i=1,iip1
	2427	ps(i,j) = ps(i,j) *
	2428	. exp((phisold(i,j)-phis(i,j))/(temp * r))
	2429	enddo
	2430	enddo
	2431	c periodicity of surface ps in longitude
	2432	do j=1,jjp1
	2433	do i=1,iip1
	2434	ps(i,j) = ps(i,j) *
	2435	. exp((phisold(i,j)-phis(i,j))/(temp * r))
	2436	enddo
	2437	enddo
	2438	c periodicity of surface ps in longitude
	2439	do j=1,jjp1
	2440	ps(1,j) = ps(iip1,j)
	2441	end do
	2442	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
	2443
	2444	c No flat topo
	2445	c -----------------------------------------------------
	2446	else if (modif(1:len_trim(modif)) .eq. 'toponoise') then
	2447
	2448	write(,) 'choice latitudes:'
	2449	587 read(,,iostat=ierr) val1,val2
	2450	if(ierr.ne.0) goto 587
	2451	write(,) 'choice longitudes:'
	2452	588 read(,,iostat=ierr) val3,val4
	2453	if(ierr.ne.0) goto 588
	2454	write(,) 'choice phi min max'
	2455	589 read(,,iostat=ierr) val5,val6
	2456	if(ierr.ne.0) goto 589
	2457	write(,) 'choice amplitude min max new phi'
	2458	585 read(,,iostat=ierr) val7,val8
	2459	if(ierr.ne.0) goto 585
	2460	c write(,) 'choice nb of random cases'
	2461	c586 read(,,iostat=ierr) val7
	2462	c if(ierr.ne.0) goto 586
	2463
	2464	do ig=1,ngridmx
	2465	IF ( (latfi(ig)*180./pi.ge.val1) .and.
	2466	& (latfi(ig)*180./pi.le.val2) .and.
	2467	& (lonfi(ig)*180./pi.ge.val3) .and.
	2468	& (lonfi(ig)*180./pi.le.val4) ) then
	2469
	2470	IF ((phisfi(ig).le.val6).and.
	2471	& (phisfi(ig).ge.val5)) then
	2472	CALL RANDOM_NUMBER(val9)
	2473	phisfi(ig)=val7+(val8-val7)*val9
	2474	ENDIF
	2475
	2476	ENDIF
	2477	enddo
	2478	phisold = phis
	2479	c update new phis and ps
	2480	CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
	2481	r = 8.314511E+0 *1000.E+0/mugaz
	2482	temp=37.
	2483	do j=1,jjp1
	2484	do i=1,iip1
	2485	ps(i,j) = ps(i,j) *
	2486	. exp((phisold(i,j)-phis(i,j))/(temp * r))
	2487	enddo
	2488	enddo
	2489	c periodicity of surface ps in longitude
	2490	do j=1,jjp1
	2491	do i=1,iip1
	2492	ps(i,j) = ps(i,j) *
	2493	. exp((phisold(i,j)-phis(i,j))/(temp * r))
	2494	enddo
	2495	enddo
	2496	c periodicity of surface ps in longitude
	2497	do j=1,jjp1
	2498	ps(1,j) = ps(iip1,j)
	2499	end do
	2500	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
	2501
	2502	c fill_sp : fill SP with N2 ice and adjust phisfi (flat SP)
	2503	c -----------------------------------------------------
	2504	else if (modif(1:len_trim(modif)) .eq. 'fill_sp') then
	2505
	2506	! def SP:
	2507	val2=-33. !lat1
	2508	val3=50. !lat2
	2509	val4=140. ! lon1
	2510	val5=-155. ! lon2
	2511	write(,) 'def SP :'
	2512	write(,) 'lat :',val2,val3
	2513	write(,) 'lon :',val4,'180 / -180',val5
	2514	write(,) 'choice phisfi ref SP (ex: -800):'
	2515	607 read(,,iostat=ierr) val6
	2516	if(ierr.ne.0) goto 607
	2517
	2518	write(,) 're scale the pressure :'
	2519	r = 8.314511E+0 *1000.E+0/mugaz
	2520	temp=40.
	2521	write(,) 'r, sale height temperature =',r,temp,g
	2522	qsurf=0.
	2523
	2524	write(,) 'latfi=',latfi
	2525	!CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
	2526	write(,) 'phis=',phis
	2527	write(,) 'phisfi=',phisfi
	2528	!CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
	2529	phisold = phis
	2530	write(,) 'phisold=',phisold
	2531	do ig=1,ngridmx
	2532
	2533	write(,) 'lat=',latfi(ig)*180./pi
	2534	write(,) 'lon=',lonfi(ig)*180./pi
	2535	write(,) 'phisfi=',phisfi(ig)
	2536	IF ( (latfi(ig)*180./pi.ge.val2) .and.
	2537	& (latfi(ig)*180./pi.le.val3) .and.
	2538	& ( ((lonfi(ig)*180./pi.ge.-180.) .and.
	2539	& (lonfi(ig)*180./pi.le.val5)) .or.
	2540	& ((lonfi(ig)*180./pi.ge.val4) .and.
	2541	& (lonfi(ig)*180./pi.le.180.))) ) then
	2542	write(,) 'hello SP',phisfi(ig),ig
	2543	IF (phisfi(ig).le.val6) then
	2544	qsurf(ig,igcm_n2)=(val6-phisfi(ig))/g*1000.
	2545	phisfi(ig)=val6
	2546	write(,) 'changes',phisfi(ig),qsurf(ig,igcm_n2)
	2547	ENDIF
	2548	ENDIF
	2549	end do
	2550	c update new phis and ps
	2551	CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
	2552	do j=1,jjp1
	2553	do i=1,iip1
	2554	ps(i,j) = ps(i,j) *
	2555	. exp((phisold(i,j)-phis(i,j))/(temp * r))
	2556	enddo
	2557	enddo
	2558	c periodicity of surface ps in longitude
	2559	do j=1,jjp1
	2560	ps(1,j) = ps(iip1,j)
	2561	end do
	2562
	2563	c bugch4 correct bug warm ch4 due to change of resolution
	2564	c -----------------------------------------------------
	2565	else if (modif(1:len_trim(modif)) .eq. 'bugch4') then
	2566	write(,) 'Ok we are going to try to solve this bug'
	2567	write(,) 'First extract a profile of tsoil and tsurf'
	2568	write(,) 'that you want to set as reference :'
	2569	write(,) 'tsoil_ref and tsurf_ref '
	2570	write(,) 'number of points to check '
	2571	546 read(,,iostat=ierr) jref1
	2572	if(ierr.ne.0) goto 546
	2573
	2574	write(,) 'choice acceptable tsurf range for ch4: t1,t2:'
	2575	547 read(,,iostat=ierr) val1,val2
	2576	if(ierr.ne.0) goto 547
	2577
	2578	! Check tsurf at all locations
	2579	DO ig=1,ngridmx
	2580	IF (qsurf(ig,igcm_ch4_ice).gt.0.001.and.
	2581	& qsurf(ig,igcm_n2).eq.0.) then
	2582	IF ((tsurf(ig).lt.val1) .or.
	2583	& (tsurf(ig).ge.val2)) then
	2584	write(,) 'Pb tsurf point ',ig
	2585
	2586	do jref=1,jref1
	2587	IF ((ig-jref.ge.1).and.qsurf(ig,igcm_n2).eq.0..and.
	2588	& (qsurf(int(ig-jref),igcm_ch4_ice).gt.0.001).and.
	2589	& (tsurf(ig-jref).gt.val1
	2590	& .and.tsurf(ig-jref).lt.val2)) then
	2591	tsurf(ig)=tsurf(int(ig-jref))
	2592	DO l=1,nsoilmx
	2593	tsoil(ig,l)=tsoil(int(ig-jref),l)
	2594	ENDDO
	2595	goto 548
	2596	ELSEIF ((ig+jref.le.ngridmx).and.
	2597	& qsurf(ig,igcm_n2).eq.0..and.
	2598	& (qsurf(int(ig+jref),igcm_ch4_ice).gt.0.001).and.
	2599	& (tsurf(ig+jref).gt.val1
	2600	& .and.tsurf(ig+jref).lt.val2)) then
	2601	tsurf(ig)=tsurf(int(ig+jref))
	2602	DO l=1,nsoilmx
	2603	tsoil(ig,l)=tsoil(int(ig+jref),l)
	2604	ENDDO
	2605	goto 548
	2606	ENDIF
	2607	enddo
	2608	548 write(,) 'found point with ch4'
	2609	ENDIF
	2610	ENDIF
	2611	END DO
	2612
	2613	c bugn2 correct bug warm n2 due to change of resolution
	2614	c -----------------------------------------------------
	2615	else if (modif(1:len_trim(modif)) .eq. 'bugn2') then
	2616	write(,) 'Ok we are going to try to solve this bug'
	2617	write(,) 'First extract a profile of tsoil and tsurf'
	2618	write(,) 'that you want to set as reference :'
	2619	write(,) 'tsoil_ref and tsurf_ref '
	2620	write(,) 'number of points to check '
	2621	544 read(,,iostat=ierr) jref1
	2622	if(ierr.ne.0) goto 544
	2623
	2624	write(,) 'choice acceptable tsurf range for n2: t1,t2:'
	2625	540 read(,,iostat=ierr) val1,val2
	2626	if(ierr.ne.0) goto 540
	2627
	2628	! Check tsurf at all locations
	2629	DO ig=1,ngridmx
	2630	IF (qsurf(ig,1).gt.0.001) then
	2631	IF ((tsurf(ig).lt.val1) .or.
	2632	& (tsurf(ig).ge.val2)) then
	2633	write(,) 'Pb tsurf point ',ig
	2634
	2635	! IF low topo et peu/pas de N2: add ch4, CO, N2
	2636	do jref=1,jref1
	2637	IF ((ig-jref.ge.1).and.
	2638	& (qsurf(int(ig-jref),igcm_n2).gt.0.001).and.
	2639	& (tsurf(ig-jref).gt.val1
	2640	& .and.tsurf(ig-jref).lt.val2)) then
	2641	tsurf(ig)=tsurf(int(ig-jref))
	2642	DO l=1,nsoilmx
	2643	tsoil(ig,l)=tsoil(int(ig-jref),l)
	2644	ENDDO
	2645	goto 541
	2646	ELSEIF ((ig+jref.le.ngridmx).and.
	2647	& (qsurf(int(ig+jref),igcm_n2).gt.0.001).and.
	2648	& (tsurf(ig+jref).gt.val1
	2649	& .and.tsurf(ig+jref).lt.val2)) then
	2650	tsurf(ig)=tsurf(int(ig+jref))
	2651	DO l=1,nsoilmx
	2652	tsoil(ig,l)=tsoil(int(ig+jref),l)
	2653	ENDDO
	2654	goto 541
	2655	ENDIF
	2656	enddo
	2657	541 write(,) 'found point with n2'
	2658	ENDIF
	2659	ENDIF
	2660	END DO
	2661
	2662	! second tour
	2663	! DO ig=1,ngridmx
	2664	! IF (qsurf(ig,1).gt.0.001) then
	2665	! IF ((tsurf(ig).lt.val1) .or.
	2666	! & (tsurf(ig).ge.val2)) then
	2667	! ! IF low topo et peu/pas de N2: add ch4, CO, N2
	2668	! IF ((ig-1.lt.1).or.(ig+1.gt.ngridmx)) then
	2669	! write(,) 'pole : doing nothing'
	2670	! ELSEIF (qsurf(ig-1,igcm_n2).gt.0.001) then
	2671	! tsurf(ig)=tsurf(ig-1)
	2672	! DO l=1,nsoilmx
	2673	! tsoil(ig,l)=tsoil(ig-1,l)
	2674	! ENDDO
	2675	! ELSEIF (qsurf(ig+1,igcm_n2).gt.0.001) then
	2676	! tsurf(ig)=tsurf(ig+1)
	2677	! DO l=1,nsoilmx
	2678	! tsoil(ig,l)=tsoil(ig+1,l)
	2679	! ENDDO
	2680	! ENDIF
	2681	! ENDIF
	2682	! ENDIF
	2683
	2684	! END DO
	2685
	2686	c flatnosp : flat topo outside a specific terrain (SP)
	2687	c -----------------------------------------------------
	2688	else if (modif(1:len_trim(modif)) .eq. 'flatnosp') then
	2689
	2690	write(,) 'Choice of topography (m) below which we keep ?'
	2691	551 read(,,iostat=ierr) val
	2692	if(ierr.ne.0) goto 551
	2693	write(,) 'gravity g is : ',g
	2694	geop=g*val
	2695	write(,) 'Choice of minimum amount of N2 ice we keep ?'
	2696	552 read(,,iostat=ierr) val2
	2697	if(ierr.ne.0) goto 552
	2698
	2699	write(,) 'phis=',phis
	2700	write(,) 'phisfi=',phisfi
	2701	do ig=1,ngridmx
	2702
	2703	IF ( (qsurf(ig,1).le.val2) .and.
	2704	& (phisfi(ig).gt.geop) ) then
	2705	write(,) 'hello SP',phisfi(ig),ig
	2706	phisfi(ig)=0.
	2707	ENDIF
	2708	end do
	2709
	2710	phisold = phis
	2711	write(,) 're scale the pressure :'
	2712	r = 8.314511E+0 *1000.E+0/mugaz
	2713	temp=40.
	2714	write(,) 'r, sale height temperature =',r,temp,g
	2715	c update new phis and ps
	2716	CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
	2717	do j=1,jjp1
	2718	do i=1,iip1
	2719	ps(i,j) = ps(i,j) *
	2720	. exp((phisold(i,j)-phis(i,j))/(temp * r))
	2721	enddo
	2722	enddo
	2723	c periodicity of surface ps in longitude
	2724	do j=1,jjp1
	2725	ps(1,j) = ps(iip1,j)
	2726	end do
	2727
	2728	c flatregion : flat topo specific to region
	2729	c -----------------------------------------------------
	2730	else if (modif(1:len_trim(modif)) .eq. 'flatregion') then
	2731
	2732	write(,) 'Choice band : lat1 and lat2 ?'
	2733	553 read(,,iostat=ierr) val,val2
	2734	if(ierr.ne.0) goto 553
	2735	write(,) 'Choice band : lon1 and lon2 ?'
	2736	554 read(,,iostat=ierr) val3,val4
	2737	if(ierr.ne.0) goto 554
	2738	write(,) 'Choice of topography range ?'
	2739	555 read(,,iostat=ierr) val5,val6
	2740	if(ierr.ne.0) goto 555
	2741	write(,) 'Choice topo ?'
	2742	556 read(,,iostat=ierr) val7
	2743	if(ierr.ne.0) goto 556
	2744
	2745	write(,) 'phis=',phis
	2746	write(,) 'phisfi=',phisfi
	2747	do ig=1,ngridmx
	2748	IF ( (latfi(ig)*180./pi.ge.val) .and.
	2749	& (latfi(ig)*180./pi.le.val2) .and.
	2750	& (lonfi(ig)*180./pi.ge.val3) .and.
	2751	& (lonfi(ig)*180./pi.le.val4) ) then
	2752
	2753	IF ( (phisfi(ig).ge.val5) .and.
	2754	& (phisfi(ig).le.val6) ) then
	2755	write(,) 'hello topo',phisfi(ig),ig
	2756	phisfi(ig)=val7
	2757	ENDIF
	2758	ENDIF
	2759	end do
	2760
	2761	r = 8.314511E+0 *1000.E+0/mugaz
	2762	temp=40.
	2763	phisold = phis
	2764	c update new phis and ps
	2765	CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
	2766	do j=1,jjp1
	2767	do i=1,iip1
	2768	ps(i,j) = ps(i,j) *
	2769	. exp((phisold(i,j)-phis(i,j))/(temp * r))
	2770	enddo
	2771	enddo
	2772	c periodicity of surface ps in longitude
	2773	do j=1,jjp1
	2774	ps(1,j) = ps(iip1,j)
	2775	end do
	2776
	2777	c changetopo
	2778	c -----------------------------------------------------
	2779	else if (modif(1:len_trim(modif)) .eq. 'changetopo') then
	2780
	2781	write(,) 'Choice band : lat1 and lat2 ?'
	2782	573 read(,,iostat=ierr) val,val2
	2783	if(ierr.ne.0) goto 573
	2784	write(,) 'Choice band : lon1 and lon2 ?'
	2785	574 read(,,iostat=ierr) val3,val4
	2786	if(ierr.ne.0) goto 574
	2787	write(,) 'Choice of topography range ?'
	2788	575 read(,,iostat=ierr) val5,val6
	2789	if(ierr.ne.0) goto 575
	2790	write(,) 'Choice change topo factor?'
	2791	576 read(,,iostat=ierr) val7
	2792	if(ierr.ne.0) goto 576
	2793	write(,) 'Choice change topo add?'
	2794	577 read(,,iostat=ierr) val8
	2795	if(ierr.ne.0) goto 577
	2796
	2797	write(,) 'phis=',phis
	2798	write(,) 'phisfi=',phisfi
	2799	do ig=1,ngridmx
	2800	IF ( (latfi(ig)*180./pi.ge.val) .and.
	2801	& (latfi(ig)*180./pi.le.val2) .and.
	2802	& (lonfi(ig)*180./pi.ge.val3) .and.
	2803	& (lonfi(ig)*180./pi.le.val4) ) then
	2804
	2805	IF ( (phisfi(ig).ge.val5) .and.
	2806	& (phisfi(ig).le.val6) ) then
	2807	write(,) 'hello topo',phisfi(ig),ig
	2808	phisfi(ig)=phisfi(ig)*val7
	2809	phisfi(ig)=phisfi(ig)+val8
	2810	ENDIF
	2811	ENDIF
	2812	end do
	2813
	2814	r = 8.314511E+0 *1000.E+0/mugaz
	2815	temp=40.
	2816	phisold = phis
	2817	c update new phis and ps
	2818	CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
	2819	do j=1,jjp1
	2820	do i=1,iip1
	2821	ps(i,j) = ps(i,j) *
	2822	. exp((phisold(i,j)-phis(i,j))/(temp * r))
	2823	enddo
	2824	enddo
	2825	c periodicity of surface ps in longitude
	2826	do j=1,jjp1
	2827	ps(1,j) = ps(iip1,j)
	2828	end do
	2829
	2830	c corrtopo: corr topo specific bug
	2831	c -----------------------------------------------------
	2832	else if (modif(1:len_trim(modif)) .eq. 'corrtopo') then
	2833
	2834	! get min max lon
	2835	print, latfi180/pi
	2836	print, '**************************'
	2837	print, '**************************'
	2838	print, '**************************'
	2839	print, '**************************'
	2840	print, '**************************'
	2841	print, lonfi180/pi
	2842	print*, 'iip1=',iip1
	2843	do ig=2,ngridmx-1
	2844	IF (lonfi(ig)*180./pi.eq.-180.) THEN
	2845	print*, lonfi(ig),lonfi(ig+iip1-2)
	2846	phisfi(ig)=(phisfi(ig+1)+phisfi(ig+iip1))/2
	2847	ENDIF
	2848	enddo
	2849	phisold = phis
	2850	r = 8.314511E+0 *1000.E+0/mugaz
	2851	temp=40.
	2852	c update new phis and ps
	2853	CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
	2854	do j=1,jjp1
	2855	do i=1,iip1
	2856	ps(i,j) = ps(i,j) *
	2857	. exp((phisold(i,j)-phis(i,j))/(temp * r))
	2858	enddo
	2859	enddo
	2860	c periodicity of surface ps in longitude
	2861	do j=1,jjp1
	2862	ps(1,j) = ps(iip1,j)
	2863	end do
	2864
	2865	c asymtopo:
	2866	c -----------------------------------------------------
	2867	else if (modif(1:len_trim(modif)) .eq. 'asymtopo') then
	2868
	2869	! get diff altitude
	2870	write(,) 'Diff N-S altitude in km (pos = N higher) ?'
	2871	591 read(,,iostat=ierr) val
	2872	if(ierr.ne.0) goto 591
	2873	write(,) 'Coeff smoot topo (small = smoother) ?'
	2874	592 read(,,iostat=ierr) val2
	2875	if(ierr.ne.0) goto 592
	2876
	2877	do ig=1,ngridmx
	2878	phisfi(ig)=phisfi(ig)+val1000.gtanh(latfi(ig)180/pi*val2)
	2879	enddo
	2880	phisold = phis
	2881	r = 8.314511E+0 *1000.E+0/mugaz
	2882	temp=40.
	2883	c update new phis and ps
	2884	CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
	2885	do j=1,jjp1
	2886	do i=1,iip1
	2887	ps(i,j) = ps(i,j) *
	2888	. exp((phisold(i,j)-phis(i,j))/(temp * r))
	2889	enddo
	2890	enddo
	2891	c periodicity of surface ps in longitude
	2892	do j=1,jjp1
	2893	ps(1,j) = ps(iip1,j)
	2894	end do
	2895
	2896	c seticeNH : set the ices in the SP crater with NH topo
	2897	c -----------------------------------------------------
	2898	else if (modif(1:len_trim(modif)) .eq. 'seticeNH') then
	2899
	2900	open(333,file='./tsoil_180_30',form='formatted')
[3390]	2901	do i=1,24
[3382]	2902	read(333,*) levsoil(i), tempsoil(i)
	2903	enddo
	2904	close(333)
	2905	open(334,file='./tsurf_180_30',form='formatted')
	2906	read(334,*) val
	2907	close(334)
	2908
	2909	write(,) 'Tsoil and tsurf ref are:'
	2910	write(,) tempsoil
	2911	write(,) 'tsurf=',val
	2912
	2913	! def SP:
	2914	val2=-43. !lat1
	2915	val3=51. !lat2
	2916	val4=143. ! lon1
	2917	val5=-158. ! lon2
	2918	val6=-150 ! phisfi min
	2919
	2920	! Rm all N2 ice outside SP
	2921	DO ig=1,ngridmx
	2922	! IF inside SP
	2923	IF ( (latfi(ig)*180./pi.ge.val2) .and.
	2924	& (latfi(ig)*180./pi.le.val3) .and.
	2925	& ( ((lonfi(ig)*180./pi.ge.-180.) .and.
	2926	& (lonfi(ig)*180./pi.le.val5)) .or.
	2927	& ((lonfi(ig)*180./pi.ge.val4) .and.
	2928	& (lonfi(ig)*180./pi.le.180.))) ) then
	2929
	2930	! IF low topo et peu/pas de N2: add ch4, CO, N2
	2931	IF ((phisfi(ig).le.val6).and.
	2932	& (qsurf(ig,igcm_n2).le.500)) then
	2933	qsurf(ig,igcm_n2)=1000.
	2934	qsurf(ig,igcm_ch4_ice)=1000.
	2935	qsurf(ig,igcm_co_ice)=1000.
	2936	tsurf(ig)=val
	2937	DO l=1,nsoilmx
	2938	tsoil(ig,l)=tempsoil(l)
	2939	ENDDO
	2940	ENDIF
	2941
	2942	! IF high topo : rm N2
	2943	IF ( (qsurf(ig,igcm_n2).ge.20.).and.
	2944	& (phisfi(ig).ge.val6) ) then
	2945	qsurf(ig,igcm_n2)=0.
	2946
	2947	ENDIF
	2948	! Mais on garde ch4 et CO en prenant la meme quantite
	2949	! qu'une autre latitude
	2950	IF ( (qsurf(ig,igcm_ch4_ice).ge.20.) .and.
	2951	& (phisfi(ig).ge.val6) ) then
	2952	jref=int(ig/jjp1)+2
	2953	qsurf(ig,igcm_ch4_ice)=qsurf(jref,igcm_ch4_ice)
	2954	ENDIF
	2955	IF ( (qsurf(ig,igcm_co_ice).ge.20.) .and.
	2956	& (phisfi(ig).ge.val6) ) then
	2957	jref=int(ig/jjp1)+2
	2958	qsurf(ig,igcm_co_ice)=qsurf(jref,igcm_co_ice)
	2959	ENDIF
	2960
	2961	ELSE ! Outside SP
	2962
	2963	IF (qsurf(ig,igcm_n2).ge.20.) then
	2964	qsurf(ig,igcm_n2)=0.
	2965	ENDIF
	2966
	2967	IF (qsurf(ig,igcm_ch4_ice).ge.20.) then
	2968	jref=int(ig/jjp1)+2
	2969	qsurf(ig,igcm_ch4_ice)=qsurf(jref,igcm_ch4_ice)
	2970	ENDIF
	2971
	2972	IF (qsurf(ig,igcm_co_ice).ge.20.) then
	2973	jref=int(ig/jjp1)+2
	2974	qsurf(ig,igcm_co_ice)=qsurf(jref,igcm_co_ice)
	2975	ENDIF
	2976
[3371]	2977	ENDIF
[3382]	2978
	2979	END DO
	2980
	2981	c 'nomountain '
	2982	c --------------
	2983	else if (modif(1:len_trim(modif)) .eq. 'nomountain') then
	2984	do j=1,jjp1
	2985	do i=1,iip1
	2986	if (phis(i,j).gt.0.1) then
	2987	phis(i,j) = 0.
	2988	ps(i,j)=ps(iim/4,j) ! assuming no topo here
	2989	endif
	2990	enddo
	2991	enddo
	2992	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
	2993
	2994	c 'relief'
	2995	c --------------
	2996	else if (modif(1:len_trim(modif)) .eq.'relief') then
	2997	write(,) "add a circular mountain/crater"
	2998	write(,) 'Center: lat lon ?'
	2999	707 read(,,iostat=ierr) lat0, lon0
	3000	if(ierr.ne.0) goto 707
	3001	if(lon0.gt.180.) lon0=lon0-360.
	3002	lat0 = lat0 * pi/180.
	3003	lon0 = lon0 * pi/180.
	3004
	3005	DO ig=1,ngridmx
	3006	if(abs(latfi(ig)-lat0).lt.pi/float(jjm) ) then
	3007	if(abs(lonfi(ig)-lon0).lt.2*pi/float(iim) ) then
	3008	ig0 = ig
	3009	end if
	3010	end if
	3011	END DO
	3012	write(,) "Reference Point to be used:"
	3013	write(,) 'ig0',ig0
	3014	write(,) 'lat=',latfi(ig0)*180./pi
	3015	write(,) 'lon=',lonfi(ig0)*180./pi
	3016
	3017	write(,) 'radius (km), height (km) negative if crater ?'
	3018	508 read(,,iostat=ierr) radm, height
	3019	if(ierr.ne.0) goto 508
	3020
	3021	write(,) 'Sale height temperature (ex:38) ?'
	3022	509 read(,,iostat=ierr) temp
	3023	if(ierr.ne.0) goto 509
	3024
	3025	r = 8.314511E+0 *1000.E+0/mugaz
	3026	do j=1,jjp1
	3027	do i=1,iip1
	3028	dist= dist_pluto(lat0,lon0,rlatu(j),rlonv(i))
	3029	if (dist.le.radm) then
	3030	phisprev= phis(i,j)
	3031	phis(i,j)=phisprev+1000.gheight*(radm-dist)/radm
	3032	write(,) 'lat=',rlatu(j)*180./pi
	3033	write(,) 'lon=',rlonv(i)*180./pi
	3034	write(,) 'dist', dist
	3035	write(,) 'z(m)=', phis(i,j)/g
	3036	ps(i,j) = ps(i,j) *
	3037	. exp((phis(i,j))/(temp * r))
	3038	end if
	3039	end do
	3040	end do
	3041
	3042	c periodicity of surface ps in longitude
	3043	do j=1,jjp1
	3044	ps(1,j) = ps(iip1,j)
	3045	end do
	3046	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
	3047
	3048	c subsoil_n2 : choice of thermal inertia values for N2 only
	3049	c ----------------------------------------------------------------
	3050	else if (modif(1:len_trim(modif)) .eq. 'subsoil_n2') then
	3051
	3052	write(,) 'New value for subsoil thermal inertia ?'
	3053	102 read(,,iostat=ierr) ith_bb
	3054	if(ierr.ne.0) goto 102
	3055	write(,) 'thermal inertia (new value):',ith_bb
	3056
	3057	write(,)'At which depth (in m.) does the ice layer begin?'
	3058	write(,)'(currently, the deepest soil layer extends down to:'
	3059	& ,layer(1),' - ',layer(nsoilmx),')'
	3060	write(,)'write 0 for uniform value for all subsurf levels?'
	3061	ierr=1
	3062	do while (ierr.ne.0)
	3063	read(,,iostat=ierr) val2
	3064	write(,)'val2:',val2,'ierr=',ierr
	3065	if (ierr.eq.0) then ! got a value, but do a sanity check
	3066	if(val2.gt.layer(nsoilmx)) then
	3067	write(,)'Depth should be less than ',layer(nsoilmx)
	3068	ierr=1
	3069	endif
	3070	if(val2.lt.layer(1)) then
	3071	if(val2.eq.0) then
	3072	write(,)'Thermal inertia set for all subsurface layers'
	3073	ierr=0
	3074	else
	3075	write(,)'Depth should be more than ',layer(1)
	3076	ierr=1
	3077	endif
	3078	endif
	3079	endif
	3080	enddo ! of do while
	3081
	3082	! find the reference index iref the depth corresponds to
	3083	if(val2.eq.0) then
	3084	iref=1
	3085	write(,)'Level selected is first level: ',layer(iref),' m'
	3086	else
	3087	do isoil=1,nsoilmx-1
	3088	if((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1)))
	3089	& then
	3090	iref=isoil
	3091	write(,)'Level selected : ',layer(isoil),' m'
	3092	exit
	3093	endif
	3094	enddo
	3095	endif
	3096
	3097	DO ig=1,ngridmx
	3098	DO j=iref,nsoilmx
	3099	c if((qsurf(ig,igcm_ch4_ice).lt.1.).and.
	3100	c & (qsurf(ig,igcm_co_ice).lt.1.))then
	3101	if (qsurf(ig,igcm_n2).gt.0.001) then
	3102	ithfi(ig,j)=ith_bb
	3103	endif
	3104	ENDDO
	3105	ENDDO
	3106
	3107	CALL gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
	3108
	3109	c subsoil_ch4 : choice of thermal inertia values for CH4 only
	3110	c ----------------------------------------------------------------
	3111	else if (modif(1:len_trim(modif)) .eq. 'subsoil_ch4') then
	3112
	3113	write(,) 'New value for subsoil thermal inertia ?'
	3114	103 read(,,iostat=ierr) ith_bb
	3115	if(ierr.ne.0) goto 103
	3116	write(,) 'thermal inertia (new value):',ith_bb
	3117
	3118	write(,)'At which depth (in m.) does the ice layer begin?'
	3119	write(,)'(currently, the deepest soil layer extends down to:'
	3120	& ,layer(1),' - ',layer(nsoilmx),')'
	3121	write(,)'write 0 for uniform value for all subsurf levels?'
	3122	ierr=1
	3123	do while (ierr.ne.0)
	3124	read(,,iostat=ierr) val2
	3125	write(,)'val2:',val2,'ierr=',ierr
	3126	if (ierr.eq.0) then ! got a value, but do a sanity check
	3127	if(val2.gt.layer(nsoilmx)) then
	3128	write(,)'Depth should be less than ',layer(nsoilmx)
	3129	ierr=1
	3130	endif
	3131	if(val2.lt.layer(1)) then
	3132	if(val2.eq.0) then
	3133	write(,)'Thermal inertia set for all subsurface layers'
	3134	ierr=0
	3135	else
	3136	write(,)'Depth should be more than ',layer(1)
	3137	ierr=1
	3138	endif
	3139	endif
	3140	endif
	3141	enddo ! of do while
	3142
	3143	! find the reference index iref the depth corresponds to
	3144	if(val2.eq.0) then
	3145	iref=1
	3146	write(,)'Level selected is first level: ',layer(iref),' m'
	3147	else
	3148	do isoil=1,nsoilmx-1
	3149	if((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1)))
	3150	& then
	3151	iref=isoil
	3152	write(,)'Level selected : ',layer(isoil),' m'
	3153	exit
	3154	endif
	3155	enddo
	3156	endif
	3157
	3158	DO ig=1,ngridmx
	3159	DO j=iref,nsoilmx
	3160	if (qsurf(ig,igcm_ch4_ice).gt.0.001.and.
	3161	& qsurf(ig,igcm_n2).lt.0.001) then
	3162	ithfi(ig,j)=ith_bb
	3163	endif
	3164	ENDDO
	3165	ENDDO
	3166
	3167	CALL gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
	3168
	3169	c subsoil_alb : choice of thermal inertia values from albedo
	3170	c ----------------------------------------------------------------
	3171	else if (modif(1:len_trim(modif)) .eq. 'subsoil_alb') then
	3172
	3173	write(,) 'Choice range albedo for defining TI'
	3174	145 read(,,iostat=ierr) val,val2
	3175	if(ierr.ne.0) goto 145
	3176	write(,) 'New value for subsoil thermal inertia ?'
	3177	144 read(,,iostat=ierr) ith_bb
	3178	if(ierr.ne.0) goto 144
	3179	write(,) 'thermal inertia (new value):',ith_bb
	3180
	3181	write(,)'At which depth (in m.) does the ice layer begin?'
	3182	write(,)'(currently, the deepest soil layer extends down to:'
	3183	& ,layer(1),' - ',layer(nsoilmx),')'
	3184	write(,)'write 0 for uniform value for all subsurf levels?'
	3185	ierr=1
	3186	do while (ierr.ne.0)
	3187	read(,,iostat=ierr) val3
	3188	if (ierr.eq.0) then ! got a value, but do a sanity check
	3189	if(val3.gt.layer(nsoilmx)) then
	3190	write(,)'Depth should be less than ',layer(nsoilmx)
	3191	ierr=1
	3192	endif
	3193	if(val3.lt.layer(1)) then
	3194	if(val3.eq.0) then
	3195	write(,)'Thermal inertia set for all subsurface layers'
	3196	ierr=0
	3197	else
	3198	write(,)'Depth should be more than ',layer(1)
	3199	ierr=1
	3200	endif
	3201	endif
	3202	endif
	3203	enddo ! of do while
	3204
	3205	! find the reference index iref the depth corresponds to
	3206	if(val3.eq.0) then
	3207	iref=1
	3208	write(,)'Level selected is first level: ',layer(iref),' m'
	3209	else
	3210	do isoil=1,nsoilmx-1
	3211	if ((val3.gt.layer(isoil)).and.(val3.lt.layer(isoil+1)))
	3212	& then
	3213	iref=isoil+1
	3214	write(,)'Level selected : ',layer(isoil+1),' m'
	3215	exit
	3216	endif
	3217	enddo
	3218	endif
	3219
	3220	DO ig=1,ngridmx
	3221	IF ( (albfi(ig).ge.val) .and.
	3222	& (albfi(ig).le.val2) ) then
	3223	DO j=iref,nsoilmx
	3224	ithfi(ig,j)=ith_bb
	3225	ENDDO
	3226	ENDIF
	3227	ENDDO
	3228
	3229	CALL gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
	3230
	3231	c subsoil_all : choice of thermal inertia values (global)
	3232	c ----------------------------------------------------------------
	3233	else if (modif(1:len_trim(modif)) .eq. 'subsoil_all') then
	3234
	3235	write(,) 'New value for subsoil thermal inertia ?'
	3236	104 read(,,iostat=ierr) ith_bb
	3237	if(ierr.ne.0) goto 104
	3238	write(,) 'thermal inertia (new value):',ith_bb
	3239
	3240	write(,)'At which depth (in m.) does the ice layer begin?'
	3241	write(,)'(currently, the deepest soil layer extends down to:'
	3242	& ,layer(1),' - ',layer(nsoilmx),')'
	3243	write(,)'write 0 for uniform value for all subsurf levels?'
	3244	ierr=1
	3245	do while (ierr.ne.0)
	3246	read(,,iostat=ierr) val2
	3247	write(,)'val2:',val2,'ierr=',ierr
	3248	if (ierr.eq.0) then ! got a value, but do a sanity check
	3249	if(val2.gt.layer(nsoilmx)) then
	3250	write(,)'Depth should be less than ',layer(nsoilmx)
	3251	ierr=1
	3252	endif
	3253	if(val2.lt.layer(1)) then
	3254	if(val2.eq.0) then
	3255	write(,)'Thermal inertia set for all subsurface layers'
	3256	ierr=0
	3257	else
	3258	write(,)'Depth should be more than ',layer(1)
	3259	ierr=1
	3260	endif
	3261	endif
	3262	endif
	3263	enddo ! of do while
	3264
	3265	! find the reference index iref the depth corresponds to
	3266	if(val2.eq.0) then
	3267	iref=1
	3268	write(,)'Level selected is first level: ',layer(iref),' m'
	3269	else
	3270	do isoil=1,nsoilmx-1
	3271	!write(,)'isoil, ',isoil,val2
	3272	!write(,)'lay(i),lay(i+1):',layer(isoil),layer(isoil+1),' m'
	3273	if ((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1)))
	3274	& then
	3275	iref=isoil+1
	3276	write(,)'Level selected : ',layer(isoil+1),' m'
	3277	exit
	3278	endif
	3279	enddo
	3280	endif
	3281
	3282	DO ig=1,ngridmx
	3283	DO j=iref,nsoilmx
	3284	ithfi(ig,j)=ith_bb
	3285	ENDDO
	3286	ENDDO
	3287
	3288	CALL gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
	3289
	3290	c diurnal_TI : choice of thermal inertia values (global)
	3291	c ----------------------------------------------------------------
	3292	else if (modif(1:len_trim(modif)) .eq. 'diurnal_TI') then
	3293
	3294	write(,) 'New value for diurnal thermal inertia ?'
	3295	106 read(,,iostat=ierr) ith_bb
	3296	if(ierr.ne.0) goto 106
	3297	write(,) 'Diurnal thermal inertia (new value):',ith_bb
	3298
	3299	write(,)'At which depth (in m.) does the ice layer ends?'
	3300	write(,)'(currently, the soil layer 1 and nsoil are:'
	3301	& ,layer(1),' - ',layer(nsoilmx),')'
	3302	ierr=1
	3303	do while (ierr.ne.0)
	3304	read(,,iostat=ierr) val2
	3305	write(,)'val2:',val2,'ierr=',ierr
	3306	if (ierr.eq.0) then ! got a value, but do a sanity check
	3307	if(val2.gt.layer(nsoilmx)) then
	3308	write(,)'Depth should be less than ',layer(nsoilmx)
	3309	ierr=1
	3310	endif
	3311	if(val2.lt.layer(1)) then
	3312	write(,)'Depth should be more than ',layer(1)
	3313	ierr=1
	3314	endif
	3315	endif
	3316	enddo ! of do while
	3317
	3318	! find the reference index iref the depth corresponds to
	3319	do isoil=1,nsoilmx-1
	3320	!write(,)'isoil, ',isoil,val2
	3321	!write(,)'lay(i),lay(i+1):',layer(isoil),layer(isoil+1),' m'
	3322	if ((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1)))
	3323	& then
	3324	iref=isoil+1
	3325	write(,)'Level selected : ',layer(isoil+1),' m'
	3326	exit
	3327	endif
	3328	enddo
	3329
	3330	DO ig=1,ngridmx
	3331	DO j=1,iref
	3332	ithfi(ig,j)=ith_bb
	3333	ENDDO
	3334	ENDDO
	3335
	3336	CALL gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
	3337
	3338	c Choice surface temperature
	3339	c -----------------------------------------------------
	3340	else if (modif(1:len_trim(modif)) .eq. 'initsurf') then
	3341
	3342	write(,) 'New value for initial surface temperature ?'
	3343	105 read(,,iostat=ierr) tsurf_bb
	3344	if(ierr.ne.0) goto 105
	3345	write(,) 'new surface temperature (K):',tsurf_bb
	3346	DO ig=1,ngridmx
	3347	tsurf(ig)=tsurf_bb
[3371]	3348	ENDDO
[3228]	3349
[3382]	3350	c Modify surface temperature (for GCM start)
	3351	c -----------------------------------------------------
	3352	else if (modif(1:len_trim(modif)) .eq. 'modtsurf') then
	3353
	3354	write(,) 'Choice band : lat1 and lat2 ?'
	3355	455 read(,,iostat=ierr) val,val2
	3356	if(ierr.ne.0) goto 455
	3357	write(,) 'Choice band : lon1 and lon2 ?'
	3358	456 read(,,iostat=ierr) val4,val5
	3359	if(ierr.ne.0) goto 456
	3360	write(,) 'Choice topo min max '
	3361	473 read(,,iostat=ierr) val9,val10
	3362	if(ierr.ne.0) goto 473
	3363	write(,) 'Choice tsurf min max '
	3364	474 read(,,iostat=ierr) val11,val12
	3365	if(ierr.ne.0) goto 474
	3366	write(,) 'Choice multiplicative factor'
	3367	457 read(,,iostat=ierr) val3
	3368	if(ierr.ne.0) goto 457
	3369	write(,) 'Choice additional tsurf'
	3370	458 read(,,iostat=ierr) val6
	3371	if(ierr.ne.0) goto 458
	3372	write(,) 'Choice multiplicative factor tsoil'
	3373	459 read(,,iostat=ierr) val7
	3374	if(ierr.ne.0) goto 459
	3375	write(,) 'Choice additional tsoil'
	3376	469 read(,,iostat=ierr) val8
	3377	if(ierr.ne.0) goto 469
	3378
	3379	DO ig=1,ngridmx
	3380	IF ( (latfi(ig)*180./pi.ge.val) .and.
	3381	& (latfi(ig)*180./pi.le.val2) .and.
	3382	& (lonfi(ig)*180./pi.ge.val4) .and.
	3383	& (lonfi(ig)*180./pi.le.val5) .and.
	3384	& (phisfi(ig).ge.val9) .and.
	3385	& (phisfi(ig).lt.val10) .and.
	3386	& (tsurf(ig).ge.val11) .and.
	3387	& (tsurf(ig).lt.val12) ) then
	3388
	3389	tsurf(ig)=tsurf(ig)*val3
	3390	tsurf(ig)=tsurf(ig)+val6
	3391	DO l=1,nsoilmx
	3392	tsoil(ig,l)=tsoil(ig,l)*val7
	3393	tsoil(ig,l)=tsoil(ig,l)+val8
	3394	ENDDO
	3395	ENDIF
	3396	ENDDO
	3397
	3398	c copy latitudes tsurf / tsoil
	3399	c -----------------------------------------------------
	3400	else if (modif(1:len_trim(modif)) .eq. 'copylat') then
	3401
	3402	write(,) 'all latitudes : ',rlatu*180./pi
	3403	write(,) 'Choice band to be modified lat1 ?'
	3404	465 read(,,iostat=ierr) val
	3405	if(ierr.ne.0) goto 465
	3406	write(,) 'Choice band copied lat2 ?'
	3407	466 read(,,iostat=ierr) val2
	3408	if(ierr.ne.0) goto 466
	3409	write(,) 'Choice multiplicative factor'
	3410	467 read(,,iostat=ierr) val3
	3411	if(ierr.ne.0) goto 467
	3412	write(,) 'Choice additional tsurf'
	3413	468 read(,,iostat=ierr) val4
	3414	if(ierr.ne.0) goto 468
	3415
	3416	j=1
	3417	DO ig=1,ngridmx
	3418	IF ((latfi(ig)*180./pi.lt.val2+180./iip1) .and.
	3419	& (latfi(ig)*180./pi.gt.val2-180./iip1)) then
	3420	randtab(j)=ig
	3421	j=j+1
	3422	ENDIF
	3423	ENDDO
	3424	print*,j, ' latitudes found'
	3425	k=1
	3426	DO ig=1,ngridmx
	3427	IF ((latfi(ig)*180./pi.lt.val+180./iip1) .and.
	3428	& (latfi(ig)*180./pi.gt.val-180./iip1)) then
	3429	tsurf(ig)=tsurf(randtab(k))*val3
	3430	tsurf(ig)=tsurf(ig)+val4
	3431	DO l=1,nsoilmx
	3432	tsoil(ig,l)=tsoil(randtab(k),l)*val3
	3433	tsoil(ig,l)=tsoil(ig,l)+val4
	3434	ENDDO
	3435	k=k+1
	3436	ENDIF
	3437	ENDDO
	3438	print*,k, ' latitudes copied'
	3439	IF (j.ne.k) stop
	3440
	3441	c copy longitudes
	3442	c -----------------------------------------------------
	3443	else if (modif(1:len_trim(modif)) .eq. 'copylon') then
	3444
	3445	write(,) 'all longitudes : ',rlonu*180./pi
	3446	write(,) 'Choice band to be modified lon1 ?'
	3447	475 read(,,iostat=ierr) val
	3448	if(ierr.ne.0) goto 475
	3449	write(,) 'Choice band copied lon2 ?'
	3450	476 read(,,iostat=ierr) val2
	3451	if(ierr.ne.0) goto 476
	3452
	3453	j=1
	3454	DO ig=2,ngridmx-1
	3455	IF ((lonfi(ig)*180./pi.lt.183.) .and.
	3456	& (lonfi(ig)*180./pi.gt.180.)) then
	3457	randtab(j)=ig
	3458	j=j+1
	3459	print*, 'lon = ',lonfi(ig)
	3460	ENDIF
	3461	ENDDO
	3462	print*,j, ' longitudes found'
	3463	k=1
	3464	DO ig=2,ngridmx-1
	3465	IF ((lonfi(ig)*180./pi.lt.180.) .and.
	3466	& (lonfi(ig)*180./pi.gt.179.)) then
	3467	tsurf(ig)=tsurf(randtab(k))
	3468	DO l=1,nsoilmx
	3469	tsoil(ig,l)=tsoil(randtab(k),l)
	3470	ENDDO
	3471	qsurf(ig,1)=qsurf(randtab(k),1)
	3472	phisfi(ig)=phisfi(randtab(k))
	3473	k=k+1
	3474	ENDIF
	3475	ENDDO
	3476	print*,k, ' longitudes copied'
	3477	IF (j.ne.k) stop
	3478	CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
	3479
	3480	c copy latlon
	3481	c -----------------------------------------------------
	3482	else if (modif(1:len_trim(modif)) .eq. 'copylatlon') then
	3483
	3484	write(,) 'all longitudes : ',rlonu*180./pi
	3485	write(,) 'Choice lat/lon to be copied ?'
	3486	495 read(,,iostat=ierr) val,val2
	3487	if(ierr.ne.0) goto 495
	3488	write(,) 'Choice indx lon1 lon2 for band ?'
	3489	496 read(,,iostat=ierr) val3,val4
	3490	if(ierr.ne.0) goto 496
	3491	write(,) 'Choice latitude indx range where we copy ?'
	3492	497 read(,,iostat=ierr) val5,val6
	3493	if(ierr.ne.0) goto 497
	3494
	3495	! choice coordinate
	3496	DO ig=2,ngridmx-1
	3497	IF ( (lonfi(ig)*180./pi.gt.val2-1) .and.
	3498	& (lonfi(ig)*180./pi.lt.val2+1) .and.
	3499	& (latfi(ig)*180./pi.lt.val+1) .and.
	3500	& (latfi(ig)*180./pi.gt.val-1) ) then
	3501	ig0=ig
	3502	print,'lat/lon=',latfi(ig)180./pi,lonfi(ig)*180./pi,ig0
	3503	ENDIF
	3504	ENDDO
	3505
	3506	DO ig=2,ngridmx-1
	3507	IF ( (lonfi(ig)*180./pi.lt.val4) .and.
	3508	& (lonfi(ig)*180./pi.gt.val3) .and.
	3509	& (latfi(ig)*180./pi.lt.val6) .and.
	3510	& (latfi(ig)*180./pi.gt.val5) .and.
	3511	& (qsurf(ig,1).lt.0.001) ) then
	3512	tsurf(ig)=tsurf(ig0)
	3513	print,'corrd=',latfi(ig)180./pi,lonfi(ig)*180./pi
	3514	DO l=1,nsoilmx
	3515	tsoil(ig,l)=tsoil(ig0,l)
	3516	ENDDO
	3517	ENDIF
	3518	ENDDO
	3519
	3520	c Choice surface temperature depending on N2 ice distribution
	3521	c -----------------------------------------------------
	3522	else if (modif(1:len_trim(modif)) .eq. 'tsurfice') then
	3523
	3524	write(,) 'Initial soil and surface temp below n2 ?'
	3525	107 read(,,iostat=ierr) tsurf_bb
	3526	if(ierr.ne.0) goto 107
	3527	write(,) 'new surface/soil temp N2 (K):',tsurf_bb
	3528	write(,) 'Initial soil and surface temp when no n2 ?'
	3529	108 read(,,iostat=ierr) tsurf_bb2
	3530	if(ierr.ne.0) goto 108
	3531	write(,) 'new surface/soil temp when no n2 (K):',tsurf_bb2
	3532	DO ig=1,ngridmx
	3533	if (qsurf(ig,igcm_n2).gt.0.001) then
	3534	tsurf(ig)=tsurf_bb
	3535	do l=1,nsoilmx
	3536	tsoil(ig,l) = tsurf_bb
	3537	end do
	3538	else if (tsurf_bb2.ne.0) then
	3539	tsurf(ig)=tsurf_bb2
	3540	do l=1,nsoilmx
	3541	tsoil(ig,l) = tsurf_bb2
	3542	end do
	3543	endif
	3544	ENDDO
	3545
	3546	c read an albedo map
	3547	c -----------------------------------------------------
	3548	else if (modif(1:len_trim(modif)) .eq. 'albedomap') then
	3549
	3550	! Get field 2D
	3551	fichnom = 'albedo.nc'
	3552	ierr = NF_OPEN (fichnom, NF_NOWRITE,nid_fi_input)
	3553	IF (ierr.NE.NF_NOERR) THEN
	3554	write(6,*)' Problem opening albedo file:',fichnom
	3555	write(6,*)' ierr = ', ierr
	3556	CALL ABORT
	3557	ENDIF
	3558
	3559	ierr = NF_INQ_VARID (nid_fi_input, trim("albedodat"),
	3560	& nvarid_input)
	3561	IF (ierr .NE. NF_NOERR) THEN
	3562	PRINT*, "Could not find asked variable in albedo.nc"
	3563	CALL abort
	3564	ENDIF
	3565
	3566	#ifdef NC_DOUBLE
	3567	ierr = NF_GET_VAR_DOUBLE(nid_fi_input,nvarid_input,
	3568	& field_input)
	3569	#else
	3570	ierr = NF_GET_VAR_REAL(nid_fi_input,nvarid_input,field_input)
	3571	#endif
	3572	IF (ierr .NE. NF_NOERR) THEN
	3573	PRINT*, "Could not get asked variable in albedo.nc"
	3574	CALL abort
	3575	ELSE
	3576	PRINT*, "Got variable in albedo.nc"
	3577	ENDIF
	3578
	3579	DO ig=1,ngridmx
	3580	albfi(ig)=field_input(ig)
	3581	ENDDO
	3582
	3583	c slope : add slope on all longitudes
	3584	c -----------------------------------------------------
	3585	else if (modif(1:len_trim(modif)) .eq. 'slope') then
	3586
	3587	write(,) 'choice latitude alt minimum & altitude value (m):'
	3588	603 read(,,iostat=ierr) val1,val2
	3589	if(ierr.ne.0) goto 603
	3590	write(,) 'choice latitude alt maximum & altitude value (m):'
	3591	604 read(,,iostat=ierr) val3,val4
	3592	if(ierr.ne.0) goto 604
	3593
	3594	write(,) 're scale the pressure :'
	3595	r = 8.314511E+0 *1000.E+0/mugaz
	3596	temp=40.
	3597	write(,) 'r, sale height temperature =',r,temp
	3598
	3599	do j=1,jjp1
	3600	do i=1,iip1
	3601	phisprev= phis(i,j)
	3602	IF ( ( (rlatu(j)*180./pi.ge.val1) .and.
	3603	& (rlatu(j)*180./pi.le.val3) ) .or.
	3604	& ( (rlatu(j)*180./pi.le.val1) .and.
	3605	& (rlatu(j)*180./pi.ge.val3) )) then
	3606
	3607	val5=abs(val2-val4)/abs(val1-val3)*
	3608	& abs(val1-rlatu(j)*180./pi)+val2
	3609	phis(i,j)=val5*g
	3610	ps(i,j) = ps(i,j) *
	3611	. exp((phisprev-phis(i,j))/(temp * r))
	3612	ENDIF
	3613	end do
	3614	end do
	3615	c periodicity of surface ps in longitude
	3616	do j=1,jjp1
	3617	ps(1,j) = ps(iip1,j)
	3618	end do
	3619	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
	3620
	3621	c digsp : change topography of SP
	3622	c -----------------------------------------------------
	3623	else if (modif(1:len_trim(modif)) .eq. 'digsp') then
	3624
	3625	write(,) 'choice latitudes:'
	3626	517 read(,,iostat=ierr) val1,val2
	3627	if(ierr.ne.0) goto 517
	3628	write(,) 'choice longitudes:'
	3629	518 read(,,iostat=ierr) val3,val4
	3630	if(ierr.ne.0) goto 518
	3631	write(,) 'choice phi limite'
	3632	519 read(,,iostat=ierr) val5
	3633	if(ierr.ne.0) goto 519
	3634	write(,) 'choice change alt (m)'
	3635	520 read(,,iostat=ierr) val6
	3636	if(ierr.ne.0) goto 520
	3637
	3638	write(,) 're scale the pressure :'
	3639	r = 8.314511E+0 *1000.E+0/mugaz
	3640	temp=40.
	3641	write(,) 'r, sale height temperature =',r,temp
	3642
	3643	do j=1,jjp1
	3644	do i=1,iip1
	3645	phisprev= phis(i,j)
	3646	IF ( ( (rlatu(j)*180./pi.ge.val1) .and.
	3647	& (rlatu(j)*180./pi.le.val2) ) .and.
	3648	& ( (rlonv(j)*180./pi.ge.val3) .and.
	3649	& (rlonv(j)*180./pi.le.val4) ) .and.
	3650	& (phis(i,j).le.val5)) then
	3651
	3652	phis(i,j)=phis(i,j)+val6*g
	3653	ps(i,j) = ps(i,j) *
	3654	. exp((phisprev-phis(i,j))/(temp * r))
	3655	ENDIF
	3656	end do
	3657	end do
	3658	c periodicity of surface ps in longitude
	3659	do j=1,jjp1
	3660	ps(1,j) = ps(iip1,j)
	3661	end do
	3662	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
	3663
	3664	c copy field 2D
	3665	c -----------------------------------------------------
	3666	else if (modif(1:len_trim(modif)) .eq. 'copytsoil') then
	3667
	3668	! Get field 2D
	3669	fichnom = 'startfi_input.nc'
	3670	ierr = NF_OPEN (fichnom, NF_NOWRITE,nid_fi_input)
	3671	IF (ierr.NE.NF_NOERR) THEN
	3672	write(6,*)' Problem opening file:',fichnom
	3673	write(6,*)' ierr = ', ierr
	3674	CALL ABORT
	3675	ENDIF
	3676
	3677	! Choice variable to be copied
	3678	c write(,) 'Choice variable to be copied ?'
	3679	c515 read(*,fmt='(a20)',iostat=ierr) modif
	3680	c if(ierr.ne.0) goto 515
	3681	c write(,) 'variable asked :',trim(modif)
	3682
	3683	ierr = NF_INQ_VARID (nid_fi_input, trim("tsurf"),nvarid_input)
	3684	IF (ierr .NE. NF_NOERR) THEN
	3685	PRINT*, "Could not find asked variable in startfi_input.nc"
	3686	CALL abort
	3687	ENDIF
	3688	ierr = NF_INQ_VARID (nid_fi_input, trim("tsoil"),
	3689	& nvarid_inputs)
	3690	IF (ierr .NE. NF_NOERR) THEN
	3691	PRINT*, "Could not find asked variable s in startfi_input.nc"
	3692	CALL abort
	3693	ENDIF
	3694
	3695	#ifdef NC_DOUBLE
	3696	ierr = NF_GET_VAR_DOUBLE(nid_fi_input,nvarid_input,
	3697	& field_input)
	3698	#else
	3699	ierr = NF_GET_VAR_REAL(nid_fi_input,nvarid_input,field_input)
	3700	#endif
	3701	IF (ierr .NE. NF_NOERR) THEN
	3702	PRINT*, "Could not get asked variable in startfi_input.nc"
	3703	CALL abort
	3704	ELSE
	3705	PRINT*, "Got variable in startfi_input.nc"
	3706	ENDIF
	3707	#ifdef NC_DOUBLE
	3708	ierr=NF_GET_VAR_DOUBLE(nid_fi_input,nvarid_inputs,
	3709	& field_inputs)
	3710	#else
	3711	ierr=NF_GET_VAR_REAL(nid_fi_input,nvarid_inputs,field_inputs)
	3712	#endif
	3713	IF (ierr .NE. NF_NOERR) THEN
	3714	PRINT*, "Could not get asked variable s in startfi_input.nc"
	3715	CALL abort
	3716	ELSE
	3717	PRINT*, "Got variable s in startfi_input.nc"
	3718	ENDIF
	3719
	3720	write(,) 'Choice domain lon1 lon2 ?'
	3721	525 read(,,iostat=ierr) val,val2
	3722	if(ierr.ne.0) goto 525
	3723	write(,) 'Choice domain lat1 lat2 ?'
	3724	526 read(,,iostat=ierr) val3,val4
	3725	if(ierr.ne.0) goto 526
	3726	write(,) 'No change if N2 ice (0) / Change (1) ?'
	3727	527 read(,,iostat=ierr) val5
	3728	if(ierr.ne.0) goto 527
	3729
	3730	! Loop
	3731	DO ig=1,ngridmx
	3732	IF ( (lonfi(ig)*180./pi.ge.val) .and.
	3733	& (lonfi(ig)*180./pi.le.val2) .and.
	3734	& (latfi(ig)*180./pi.ge.val3) .and.
	3735	& (latfi(ig)*180./pi.le.val4) ) then
	3736	if (qsurf(ig,igcm_n2).lt.0.001.or.val5.eq.1) then
	3737	tsurf(ig)=field_input(ig)
	3738	do l=1,nsoilmx
[3390]	3739	!tsoil(ig,l) = field_input(ig)
	3740	tsoil(ig,l) = field_inputs(ig,l)
[3382]	3741	end do
	3742	endif
	3743	ENDIF
	3744	ENDDO
	3745
	3746	c modify ice distribution depending on albedo
	3747	c -----------------------------------------------------
	3748	else if (modif(1:len_trim(modif)) .eq. 'mod_distrib_ice') then
	3749
	3750	write(,) 'Choice range albedo for CH4 ice'
	3751	220 read(,,iostat=ierr) val,val2
	3752	if(ierr.ne.0) goto 220
	3753	write(,) 'Choice range albedo for N2 ice'
	3754	221 read(,,iostat=ierr) val3,val4
	3755	if(ierr.ne.0) goto 221
	3756	write(,) 'Choice extra mass of CH4 ice'
	3757	222 read(,,iostat=ierr) val5
	3758	if(ierr.ne.0) goto 222
	3759	write(,) 'Choice extra mass of N2 ice'
	3760	223 read(,,iostat=ierr) val6
	3761	if(ierr.ne.0) goto 223
	3762
	3763	DO ig=1,ngridmx
	3764	IF ( (albfi(ig).ge.val) .and.
	3765	& (albfi(ig).le.val2) ) then
	3766	qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)+val5
	3767	ENDIF
	3768	IF ( (albfi(ig).ge.val3) .and.
	3769	& (albfi(ig).le.val4) ) then
	3770	qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)+val6
	3771	ENDIF
	3772	ENDDO
	3773
[3371]	3774	c ----------------------------------------------------------
	3775	c 'reservoir' : increase or decrase ices reservoir by factor
	3776	c ----------------------------------------------------------
	3777	else if (modif(1:len_trim(modif)).eq.'reservoir') then
	3778	write(,) "Increase/Decrease N2 reservoir by factor :"
	3779	read(,) val
	3780	write(,) "Increase/Decrease CH4 reservoir by factor :"
	3781	read(,) val2
	3782	write(,) "Increase/Decrease CO reservoir by factor :"
	3783	read(,) val3
[3228]	3784
[3371]	3785	DO ig=1,ngridmx
	3786	qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)*val
	3787	qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)*val2
	3788	qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)*val3
	3789	ENDDO
[3228]	3790
[3371]	3791	c --------------------------------------------------------
	3792	c 'plutomap' : initialize pluto ices distribution
	3793	c --------------------------------------------------------
	3794	else if (modif(1:len_trim(modif)).eq.'plutomap') then
[3228]	3795
[3371]	3796	write(,) 'pluto_map.dat has to be in your simulation file !!'
	3797	open(27,file='pluto_map.dat')
	3798	N2_pluto_dat(:,:) =0.
	3799	CH4_pluto_dat(:,:) =0.
	3800	CO_pluto_dat(:,:) =0.
[3228]	3801
[3371]	3802	! Dimension file pluto_map
	3803	IF (jm_plu.ne.180) then
	3804	write(,) 'Newstart : you should set jm_plu to 180'
	3805	stop
	3806	ENDIF
	3807	! Read values
	3808	do j=1,jm_plu+1
	3809	read(27,*) (map_pluto_dat(i,j),i=1,im_plu)
	3810	do i=1,im_plu
[3228]	3811
[3371]	3812	!!!!!! BTD_v2
	3813	if (map_pluto_dat(i,j).eq.3) then
	3814	N2_pluto_dat(i,j)=100000.
[3382]	3815	else if (map_pluto_dat(i,j).eq.4) then
[3371]	3816	CH4_pluto_dat(i,j)=100000.
[3382]	3817	else if (map_pluto_dat(i,j).eq.0) then
[3371]	3818	alb_dat(i,j)=0.3
[3382]	3819	else if (map_pluto_dat(i,j).eq.6) then
[3371]	3820	alb_dat(i,j)=0.6
[3382]	3821	else if (map_pluto_dat(i,j).eq.7) then
[3371]	3822	alb_dat(i,j)=0.1
	3823	endif
[3228]	3824
[3371]	3825	!!!!!! BTD
	3826	!if (map_pluto_dat(i,j).eq.3) then
	3827	! CH4_pluto_dat(i,j)=100000.
	3828	!endif
	3829	end do
	3830	end do
	3831	close (27)
	3832	! Interpolate
[3228]	3833
[3371]	3834	!! latitude and longitude in REindexed pluto map :
	3835	latv_plu(1)=90. *pi/180.
	3836	do j=2,jm_plu
	3837	latv_plu(j)= (90-(j-1 -0.5)180./(jm_plu-1))pi/180.
	3838	end do
	3839	latv_plu(jm_plu+1)= -90. *pi/180.
[3228]	3840
[3371]	3841	do i=1,im_plu+1
	3842	lonu_plu(i)=(-180+ (i-1)360./(im_plu) ) pi/180.
	3843	enddo
[3228]	3844
[3371]	3845	!Reindexation to shift the longitude grid like a LMD GCM grid...
	3846	do j=1,jm_plu+1
	3847	N2_pluto_rein(1,j)=(N2_pluto_dat(1,j)+
	3848	& N2_pluto_dat(im_plu,j))/2
	3849	CH4_pluto_rein(1,j)=(CH4_pluto_dat(1,j)+
	3850	& CH4_pluto_dat(im_plu,j))/2
	3851	CO_pluto_rein(1,j)=(CO_pluto_dat(1,j)+
	3852	& CO_pluto_dat(im_plu,j))/2
	3853	alb_rein(1,j)=(alb_dat(1,j)+
	3854	& alb_dat(im_plu,j))/2
	3855	do i=2,im_plu
	3856	N2_pluto_rein(i,j)=(N2_pluto_dat(i-1,j)+
	3857	& N2_pluto_dat(i,j))/2
	3858	CH4_pluto_rein(i,j)=(CH4_pluto_dat(i-1,j)+
	3859	& CH4_pluto_dat(i,j))/2
	3860	CO_pluto_rein(i,j)=(CO_pluto_dat(i-1,j)+
	3861	& CO_pluto_dat(i,j))/2
	3862	alb_rein(i,j)=(alb_dat(i-1,j)+
	3863	& alb_dat(i,j))/2
	3864	end do
	3865	N2_pluto_rein(im_plu+1,j) = N2_pluto_rein(1,j)
	3866	CH4_pluto_rein(im_plu+1,j) = CH4_pluto_rein(1,j)
	3867	CO_pluto_rein(im_plu+1,j) = CO_pluto_rein(1,j)
	3868	alb_rein(im_plu+1,j) = alb_rein(1,j)
	3869	end do
[3228]	3870
[3371]	3871	call interp_horiz(N2_pluto_rein,N2_pluto_gcm,
	3872	& im_plu,jm_plu,iim,jjm,1,lonu_plu,latv_plu,rlonu,rlatv)
	3873	call interp_horiz(CH4_pluto_rein,CH4_pluto_gcm,
	3874	& im_plu,jm_plu,iim,jjm,1,lonu_plu,latv_plu,rlonu,rlatv)
	3875	call interp_horiz(CO_pluto_rein,CO_pluto_gcm,
	3876	& im_plu,jm_plu,iim,jjm,1,lonu_plu,latv_plu,rlonu,rlatv)
	3877	call interp_horiz(alb_rein,alb_gcm,
	3878	& im_plu,jm_plu,iim,jjm,1,lonu_plu,latv_plu,rlonu,rlatv)
[3228]	3879
[3371]	3880	! Switch dyn grid to fi grid
	3881	qsurf_tmp(:,:) =0.
	3882	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,
	3883	& N2_pluto_gcm,qsurf_tmp(1,igcm_n2))
	3884	if (igcm_ch4_ice.ne.0) then
	3885	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,
	3886	& CH4_pluto_gcm,qsurf_tmp(1,igcm_ch4_ice))
	3887	endif
	3888	if (igcm_co_ice.ne.0) then
	3889	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,
	3890	& CO_pluto_gcm,qsurf_tmp(1,igcm_co_ice))
	3891	endif
	3892	alb=alb_gcm
	3893	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb_gcm,albfi)
	3894	!print*, 'N2dat=',N2_pluto_gcm
	3895	!print*, 'COdat=',CO_pluto_gcm
	3896	!print*, 'CH4dat=',CH4_pluto_gcm
	3897	print*, 'N2dat=',qsurf_tmp(:,igcm_n2)
	3898	print*, 'COdat=',qsurf_tmp(:,igcm_co_ice)
	3899	print*, 'CH4dat=',qsurf_tmp(:,igcm_ch4_ice)
	3900
	3901	! Fill
	3902	DO ig=1,ngridmx
	3903	qsurf(ig,igcm_n2)=qsurf(ig,igcm_n2)+qsurf_tmp(ig,igcm_n2)
	3904	if (igcm_ch4_ice.ne.0) then
	3905	qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)+
	3906	& qsurf_tmp(ig,igcm_ch4_ice)
	3907	endif
	3908	if (igcm_co_ice.ne.0) then
	3909	qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)+
	3910	& qsurf_tmp(ig,igcm_co_ice)
	3911	endif
	3912	ENDDO
	3913
[3184]	3914	endif
[3198]	3915
[3184]	3916	enddo ! of do ! infinite loop on liste of changes
	3917
	3918	999 continue
	3919
	3920	c=======================================================================
	3921	c Correct pressure on the new grid (menu 0)
	3922	c=======================================================================
	3923
	3924	if ((choix_1.eq.0).and.(.not.flagps0)) then
	3925	r = 1000.*8.31/mugaz
	3926
	3927	do j=1,jjp1
	3928	do i=1,iip1
	3929	ps(i,j) = ps(i,j) *
	3930	. exp((phisold_newgrid(i,j)-phis(i,j)) /
	3931	. (t(i,j,1) * r))
	3932	end do
	3933	end do
	3934
	3935	c periodicite de ps en longitude
	3936	do j=1,jjp1
	3937	ps(1,j) = ps(iip1,j)
	3938	end do
	3939	end if
	3940
	3941	c=======================================================================
	3942	c Initialisation de la physique / ecriture de newstartfi :
	3943	c=======================================================================
	3944
	3945	CALL inifilr
	3946	CALL pression(ip1jmp1, ap, bp, ps, p3d)
	3947
	3948	c-----------------------------------------------------------------------
	3949	c Initialisation pks:
	3950	c-----------------------------------------------------------------------
	3951
	3952	CALL exner_hyb(ip1jmp1, ps, p3d, pks, pk, pkf)
	3953	! Calcul de la temperature potentielle teta
	3954
	3955	if (flagtset) then
	3956	DO l=1,llm
	3957	DO j=1,jjp1
	3958	DO i=1,iim
	3959	teta(i,j,l) = Tset(i,j,l) * cpp/pk(i,j,l)
	3960	ENDDO
	3961	teta (iip1,j,l)= teta (1,j,l)
	3962	ENDDO
	3963	ENDDO
	3964	else if (choix_1.eq.0) then
	3965	DO l=1,llm
	3966	DO j=1,jjp1
	3967	DO i=1,iim
	3968	teta(i,j,l) = t(i,j,l) * cpp/pk(i,j,l)
	3969	ENDDO
	3970	teta (iip1,j,l)= teta (1,j,l)
	3971	ENDDO
	3972	ENDDO
	3973	endif
	3974
	3975	C Calcul intermediaire
	3976	c
	3977	if (choix_1.eq.0) then
	3978	CALL massdair( p3d, masse )
	3979	c
	3980	print *,' ALPHAX ',alphax
	3981
	3982	DO l = 1, llm
	3983	DO i = 1, iim
	3984	xppn(i) = aire( i, 1 ) * masse( i , 1 , l )
	3985	xpps(i) = aire( i,jjp1 ) * masse( i , jjp1 , l )
	3986	ENDDO
	3987	xpn = SUM(xppn)/apoln
	3988	xps = SUM(xpps)/apols
	3989	DO i = 1, iip1
	3990	masse( i , 1 , l ) = xpn
	3991	masse( i , jjp1 , l ) = xps
	3992	ENDDO
	3993	ENDDO
	3994	endif
	3995	phis(iip1,:) = phis(1,:)
	3996
	3997	itau=0
	3998	if (choix_1.eq.0) then
	3999	day_ini=int(date)
	4000	endif
	4001	c
	4002	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
	4003
	4004	CALL caldyn0( itau,ucov,vcov,teta,ps,masse,pk,phis ,
	4005	* phi,w, pbaru,pbarv,day_ini+time )
	4006
	4007
	4008	CALL dynredem0("restart.nc",day_ini,phis)
	4009	CALL dynredem1("restart.nc",0.0,vcov,ucov,teta,q,masse,ps)
	4010	C
	4011	C Ecriture etat initial physique
	4012	C
	4013
	4014	call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngridmx,
	4015	& llm,
	4016	& nqtot,dtphys,real(day_ini),0.0,
[3539]	4017	& cell_area,albfi,zmea,zstd,zsig,zgam,zthe)
[3184]	4018	call physdem1("restartfi.nc",nsoilmx,ngridmx,llm,nqtot,
	4019	& dtphys,real(day_ini),
[3539]	4020	& tsurf,tsoil,ithfi,emis,q2,qsurf)
[3198]	4021	! & cloudfrac,totalfrac,hice,
	4022	! & rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
[3184]	4023
	4024	c=======================================================================
	4025	c Formats
	4026	c=======================================================================
	4027
	4028	1 FORMAT(//10x,'la valeur de im =',i4,2x,'lue sur le fichier de dema
	4029	*rrage est differente de la valeur parametree iim =',i4//)
	4030	2 FORMAT(//10x,'la valeur de jm =',i4,2x,'lue sur le fichier de dema
	4031	*rrage est differente de la valeur parametree jjm =',i4//)
	4032	3 FORMAT(//10x,'la valeur de lllm =',i4,2x,'lue sur le fichier demar
	4033	*rage est differente de la valeur parametree llm =',i4//)
	4034
	4035	write(,) "newstart: All is well that ends well."
	4036
	4037	end
	4038
[3371]	4039	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	4040	function dist_pluto(lat1,lon1,lat2,lon2)
	4041	implicit none
	4042	real dist_pluto
	4043	real lat1,lon1,lat2,lon2
	4044	real dlon, dlat
	4045	real a,c
[3382]	4046	real radp
	4047	radp=1190. ! km
[3371]	4048
	4049	dlon = lon2 - lon1
	4050	dlat = lat2 - lat1
	4051	a = (sin(dlat/2))*2 + cos(lat1) cos(lat2) * (sin(dlon/2))**2
	4052	c = 2 * atan2( sqrt(a), sqrt(1-a) )
[3382]	4053	dist_pluto = radp * c
[3371]	4054	return
	4055	end
	4056

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

Download in other formats:

Original Format