Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

newstart.F @ 3603

Last change on this file since 3603 was 3175, checked in by emillour, 17 months ago
Pluto PCM: Add the old Pluto LMDZ for reference (required prior step to making an LMDZ.PLUTO using the same framework as the other physics packages). TB+EM
File size: 151.8 KB

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

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

Download in other formats:

Original Format