Context Navigation

newstart.F @ 3422

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

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

Download in other formats:

Original Format