Context Navigation

newstart.F @ 3546

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

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

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

Download in other formats:

Original Format