Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

newstart.F @ 3436

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

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

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

Download in other formats:

Original Format