Context Navigation

newstart.F @ 4025

Last change on this file since 4025 was 4025, checked in by tbertrand, 8 days ago
PLUTO PCM: Fix in newstart to allow topography to remain unchanged when using startarchive.nc TB
File size: 163.4 KB

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

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

Download in other formats:

Original Format