Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

newstart.F_ini @ 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: 58.0 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 Derniere modif : 12/03
9	c
10	c
11	c Objet: Create or modify the initial state for the LMD Mars GCM
12	c ----- (fichiers NetCDF start et startfi)
13	c
14	c
15	c=======================================================================
16
17	implicit none
18
19	#include "dimensions.h"
20	#include "dimphys.h"
21	#include "surfdat.h"
22	#include "comsoil.h"
23	#include "planete.h"
24	#include "paramet.h"
25	#include "comconst.h"
26	#include "comvert.h"
27	#include "comgeom2.h"
28	#include "control.h"
29	#include "logic.h"
30	#include "description.h"
31	#include "ener.h"
32	#include "temps.h"
33	#include "lmdstd.h"
34	#include "comdissnew.h"
35	#include "clesph0.h"
36	#include "serre.h"
37	#include "netcdf.inc"
38	#include "advtrac.h"
39	#include "tracer.h"
40	c=======================================================================
41	c Declarations
42	c=======================================================================
43
44	c Variables dimension du fichier "start_archive"
45	c------------------------------------
46	CHARACTER relief*3
47
48
49	c Variables pour les lectures NetCDF des fichiers "start_archive"
50	c--------------------------------------------------
51	INTEGER nid_dyn, nid_fi,nid,nvarid
52	INTEGER length
53	parameter (length = 100)
54	INTEGER tab0
55	INTEGER NB_ETATMAX
56	parameter (NB_ETATMAX = 100)
57
58	REAL date
59	REAL p_rad,p_omeg,p_g,p_cpp,p_mugaz,p_daysec
60
61	c Variable histoire
62	c------------------
63	REAL vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm) ! vents covariants
64	REAL phis(iip1,jjp1)
65	REAL q(iip1,jjp1,llm,nqmx) ! champs advectes
66
67	c autre variables dynamique nouvelle grille
68	c------------------------------------------
69	REAL pks(iip1,jjp1)
70	REAL w(iip1,jjp1,llm+1)
71	REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
72	! REAL dv(ip1jm,llm),du(ip1jmp1,llm)
73	! REAL dh(ip1jmp1,llm),dp(ip1jmp1)
74	REAL phi(iip1,jjp1,llm)
75
76	integer klatdat,klongdat
77	PARAMETER (klatdat=180,klongdat=360)
78
79	c Physique sur grille scalaire
80	c----------------------------
81	real zmeaS(iip1,jjp1),zstdS(iip1,jjp1)
82	real zsigS(iip1,jjp1),zgamS(iip1,jjp1),ztheS(iip1,jjp1)
83
84	c variable physique
85	c------------------
86	REAL tsurf(ngridmx) ! surface temperature
87	REAL tsoil(ngridmx,nsoilmx) ! soil temperature
88	REAL lati(ngridmx)
89	! REAL co2ice(ngridmx) ! CO2 ice layer
90	REAL emis(ngridmx) ! surface emissivity
91	real emisread ! added by RW
92	REAL qsurf(ngridmx,nqmx)
93	REAL q2(ngridmx,nlayermx+1)
94	! REAL rnaturfi(ngridmx)
95	real alb(iip1,jjp1),albfi(ngridmx) ! albedos
96	real ith(iip1,jjp1,nsoilmx),ithfi(ngridmx,nsoilmx) ! thermal inertia (3D)
97	real surfith(iip1,jjp1),surfithfi(ngridmx) ! surface thermal inertia (2D)
98	REAL latfi(ngridmx),lonfi(ngridmx),airefi(ngridmx)
99
100	INTEGER i,j,l,isoil,ig,idum
101	! real mugaz ! molar mass of the atmosphere.
102	!TB15 : in comconst.h
103
104	integer ierr
105
106	c Variables on the new grid along scalar points
107	c------------------------------------------------------
108	! REAL p(iip1,jjp1)
109	REAL t(iip1,jjp1,llm)
110	REAL tset(iip1,jjp1,llm)
111	real phisold_newgrid(iip1,jjp1)
112	! REAL tanh(ip1jmp1)
113	REAL :: teta(iip1, jjp1, llm)
114	REAL :: pk(iip1,jjp1,llm)
115	REAL :: pkf(iip1,jjp1,llm)
116	REAL :: ps(iip1, jjp1)
117	REAL :: ppa(iip1,jjp1,llm)
118	REAL :: masse(iip1,jjp1,llm)
119	REAL :: xpn,xps,xppn(iim),xpps(iim)
120	REAL :: p3d(iip1, jjp1, llm+1)
121	REAL :: beta(iip1,jjp1,llm)
122	! REAL dteta(ip1jmp1,llm)
123
124	c Variable de l'ancienne grille
125	c------------------------------
126	real time
127	real tab_cntrl(100)
128	real tab_cntrl_bis(100)
129
130	c variables diverses
131	c-------------------
132	real choix_1,choice,pp
133	character*80 fichnom
134	integer Lmodif,iq,thermo
135	character modif*20
136	real z_reel(iip1,jjp1)
137	! real z_reel(ngridmx)
138	real tsud,albsud,alb_bb,ith_bb,Tiso,Tabove,tsurf_bb,geop
139	real ptoto,pcap,patm,airetot,ptotn,patmn,psea
140	character(len=50) :: surfacefile ! "surface.nc" (or equivalent file)
141	! real ssum
142	character*1 yes
143	logical :: flagtset=.false. , flagps0=.false.
144	real val, val1, val2, val3, val4 ! to store temporary variables
145	real :: iceith=2000 ! thermal inertia of subterranean ice
146	integer iref,jref,iref1,jref1,iref2,jref2
147
148	INTEGER :: itau
149
150	INTEGER :: nq,numvanle
151	character(len=20) :: txt ! to store some text
152	integer :: count
153
154	! MONS data:
155	! real :: MONS_Hdn(iip1,jjp1) ! Hdn: %WEH=Mass fraction of H2O
156	! real :: MONS_d21(iip1,jjp1) ! ice table "depth" (in kg/m2)
157	! ! coefficient to apply to convert d21 to 'true' depth (m)
158	! real :: MONS_coeff
159	! real :: MONS_coeffS ! coeff for southern hemisphere
160	! real :: MONS_coeffN ! coeff for northern hemisphere
161	!! real,parameter :: icedepthmin=1.e-3 ! Ice begins at most at that depth
162
163	c Special Pluto Map from Lellouch & Grundy
164	c------------------------------------------
165	integer,parameter :: im_plu=59
166	integer,parameter :: jm_plu=30
167	real latv_plu(jm_plu),lonu_plu(im_plu+1)
168	real map_pluto_dat(im_plu,jm_plu+1)
169
170	real alb_pluto_dat(im_plu,jm_plu+1)
171	real N2_pluto_dat(im_plu,jm_plu+1)
172	real CH4_pluto_dat(im_plu,jm_plu+1)
173	real CO_pluto_dat(im_plu,jm_plu+1)
174	real ith_pluto_dat(im_plu,jm_plu+1)
175
176	real alb_pluto_rein(im_plu+1,jm_plu+1)
177	real N2_pluto_rein(im_plu+1,jm_plu+1)
178	real CH4_pluto_rein(im_plu+1,jm_plu+1)
179	real CO_pluto_rein(im_plu+1,jm_plu+1)
180	real ith_pluto_rein(im_plu+1,jm_plu+1)
181	real N2_pluto_gcm(iip1,jjp1)
182	real CH4_pluto_gcm(iip1,jjp1)
183	real CO_pluto_gcm(iip1,jjp1)
184	real alb_pluto_gcm(iip1,jjp1),alb_pluto_fi(ngridmx)
185	real ith_pluto_gcm(iip1,jjp1),ithf(ngridmx)
186
187	c sortie visu pour les champs dynamiques
188	c---------------------------------------
189	! INTEGER :: visuid
190	! real :: time_step,t_ops,t_wrt
191	! CHARACTER*80 :: visu_file
192
193	! preff = 2. ! Pluto ! 610. ! for Mars, instead of 101325. (Earth)
194	! pa = 0.5 ! Pluto ! 20 ! for Mars, instead of 500 (Earth)
195
196	c=======================================================================
197	c Choice of the start file(s) to use
198	c=======================================================================
199
200	write(,) 'From which kind of files do you want to create new',
201	. 'start and startfi files'
202	write(,) ' 0 - from a file start_archive'
203	write(,) ' 1 - from files start and startfi'
204
205	c-----------------------------------------------------------------------
206	c Open file(s) to modify (start or start_archive)
207	c-----------------------------------------------------------------------
208
209	DO
210	read(,,iostat=ierr) choix_1
211	if ((choix_1 /= 0).OR.(choix_1 /=1)) EXIT
212	ENDDO
213
214	c Open start_archive
215	c ~~~~~~~~~~~~~~~~~~~~~~~~~~
216	if (choix_1.eq.0) then
217
218	write(,) 'Creating start files from:'
219	write(,) './start_archive.nc'
220	write(,)
221	fichnom = 'start_archive.nc'
222	ierr = NF_OPEN (fichnom, NF_NOWRITE,nid)
223	IF (ierr.NE.NF_NOERR) THEN
224	write(6,*)' Problem opening file:',fichnom
225	write(6,*)' ierr = ', ierr
226	CALL ABORT
227	ENDIF
228	tab0 = 50
229	Lmodif = 1
230
231	c OR open start and startfi files
232	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
233	else
234	write(,) 'Creating start files from:'
235	write(,) './start.nc and ./startfi.nc'
236	write(,)
237	fichnom = 'start.nc'
238	ierr = NF_OPEN (fichnom, NF_NOWRITE,nid_dyn)
239	IF (ierr.NE.NF_NOERR) THEN
240	write(6,*)' Problem opening file:',fichnom
241	write(6,*)' ierr = ', ierr
242	CALL ABORT
243	ENDIF
244
245	fichnom = 'startfi.nc'
246	ierr = NF_OPEN (fichnom, NF_NOWRITE,nid_fi)
247	IF (ierr.NE.NF_NOERR) THEN
248	write(6,*)' Problem opening file:',fichnom
249	write(6,*)' ierr = ', ierr
250	CALL ABORT
251	ENDIF
252
253	tab0 = 0
254	Lmodif = 0
255
256	endif
257
258	c-----------------------------------------------------------------------
259	c Lecture du tableau des parametres du run (pour la dynamique)
260	c-----------------------------------------------------------------------
261
262	if (choix_1.eq.0) then
263
264	write(,) 'reading tab_cntrl START_ARCHIVE'
265	c
266	ierr = NF_INQ_VARID (nid, "controle", nvarid)
267	#ifdef NC_DOUBLE
268	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tab_cntrl)
269	#else
270	ierr = NF_GET_VAR_REAL(nid, nvarid, tab_cntrl)
271	#endif
272	c
273	else if (choix_1.eq.1) then
274
275	write(,) 'reading tab_cntrl START'
276	c
277	ierr = NF_INQ_VARID (nid_dyn, "controle", nvarid)
278	#ifdef NC_DOUBLE
279	ierr = NF_GET_VAR_DOUBLE(nid_dyn, nvarid, tab_cntrl)
280	#else
281	ierr = NF_GET_VAR_REAL(nid_dyn, nvarid, tab_cntrl)
282	#endif
283	c
284	write(,) 'reading tab_cntrl STARTFI'
285	c
286	ierr = NF_INQ_VARID (nid_fi, "controle", nvarid)
287	#ifdef NC_DOUBLE
288	ierr = NF_GET_VAR_DOUBLE(nid_fi, nvarid, tab_cntrl_bis)
289	#else
290	ierr = NF_GET_VAR_REAL(nid_fi, nvarid, tab_cntrl_bis)
291	#endif
292	c
293	do i=1,50
294	tab_cntrl(i+50)=tab_cntrl_bis(i)
295	enddo
296	endif
297
298	c TB15 :
299	c do i=1,size(tab_cntrl)
300	c write(,) i,tab_cntrl(i)
301	c enddo
302
303	c-----------------------------------------------------------------------
304	c Initialisation des constantes dynamique
305	c-----------------------------------------------------------------------
306
307	kappa = tab_cntrl(9)
308	etot0 = tab_cntrl(12)
309	ptot0 = tab_cntrl(13)
310	ztot0 = tab_cntrl(14)
311	stot0 = tab_cntrl(15)
312	ang0 = tab_cntrl(16)
313
314	write(,) "Newstart: in nc: kappa,etot0,ptot0,ztot0,stot0,ang0"
315	write(,) kappa,etot0,ptot0,ztot0,stot0,ang0
316
317	! for vertical coordinate
318	preff=tab_cntrl(18) ! reference surface pressure
319	pa=tab_cntrl(17) ! reference pressure at which coord is purely pressure
320	!NB: in start_archive files tab_cntrl(17)=tab_cntrl(18)=0
321	if (preff.eq.0) then
322	preff=2. ! Pluto
323	pa=0.5 ! Pluto
324	endif
325	write(,) "Newstart: preff=",preff," pa=",pa
326	yes=' '
327	do while ((yes.ne.'y').and.(yes.ne.'n'))
328	write(,) "Change the values of preff and pa? (y/n)"
329	read(*,fmt='(a)') yes
330	end do
331	if (yes.eq.'y') then
332	write(,)"New value of reference surface pressure preff?"
333	write(,)" (for Mars, typically preff=610)"
334	read(,) preff
335	write(,)"New value of reference pressure pa for purely"
336	write(,)"pressure levels (for hybrid coordinates)?"
337	write(,)" (for Mars, typically pa=20)"
338	read(,) pa
339	endif
340	c-----------------------------------------------------------------------
341	c Lecture du tab_cntrl et initialisation des constantes physiques
342	c - pour start: Lmodif = 0 => pas de modifications possibles
343	c (modif dans le tabfi de readfi + loin)
344	c - pour start_archive: Lmodif = 1 => modifications possibles
345	c-----------------------------------------------------------------------
346
347	if (choix_1.eq.0) then
348	call tabfi (nid,Lmodif,tab0,day_ini,lllm,p_rad,
349	. p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
350	else if (choix_1.eq.1) then
351	call tabfi (nid_fi,Lmodif,tab0,day_ini,lllm,p_rad,
352	. p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
353	endif
354
355	rad = p_rad
356	omeg = p_omeg
357	g = p_g
358	cpp = p_cpp
359	mugaz = p_mugaz
360	daysec = p_daysec
361	kappa= 8.314511E+0 1000.E+0/(mugazcpp)
362
363	c=======================================================================
364	c INITIALISATIONS DIVERSES
365	c=======================================================================
366	! Load tracer names:
367	call iniadvtrac(nq,numvanle)
368	! tnom(:) now contains tracer names
369	! Initialize global tracer indexes (stored in tracer.h)
370	call initracer()
371	! Load parameters from run.def file
372	CALL defrun_new( 99, .TRUE. )
373	CALL iniconst
374	CALL inigeom
375	idum=-1
376	idum=0
377
378	c Initialisation coordonnees /aires
379	c -------------------------------
380	! Note: rlatu(:) and rlonv(:) are commons defined in "comgeom.h"
381	! rlatu() and rlonv() are given in radians
382	latfi(1)=rlatu(1)
383	lonfi(1)=0.
384	DO j=2,jjm
385	DO i=1,iim
386	latfi((j-2)*iim+1+i)=rlatu(j)
387	lonfi((j-2)*iim+1+i)=rlonv(i)
388	ENDDO
389	ENDDO
390	latfi(ngridmx)=rlatu(jjp1)
391	lonfi(ngridmx)=0.
392	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
393
394	c=======================================================================
395	c lecture topographie, albedo, inertie thermique, relief sous-maille
396	c=======================================================================
397
398	if (choix_1.ne.1) then ! pour ne pas avoir besoin du fichier
399	! surface.dat dans le cas des start
400
401	c do while((relief(1:3).ne.'mol').AND.(relief(1:3).ne.'pla'))
402	c write(,)
403	c write(,) 'choix du relief (mola,pla)'
404	c write(,) '(Topographie MGS MOLA, plat)'
405	c read(*,fmt='(a3)') relief
406	c relief="mola"
407	c enddo
408
409	write(,) ''
410	write(,) 'Please choose:'
411	write(,) '1- use of surface file netcdf'
412	write(,) '2- define topography manually'
413	write(,) '3- skip this part'
414	read(,) choice
415
416	if (choice.eq.1) then
417
418	write(,) ' enter surface file: (e.g. "surface_mars.nc)'
419	read(*,fmt='(a50)') surfacefile
420
421	CALL datareadnc(relief,surfacefile,phis,alb,surfith,
422	& zmeaS,zstdS,zsigS,zgamS,ztheS)
423
424	else if (choice.eq.2) then
425	! specific case when not using a "surface.nc" file
426	302 write(,) 'Do you want a flat terrain? (y/n)'
427	read(*,fmt='(a)') yes
428	if (yes.eq.'y') then
429	phis(:,:)=0
430	else
431	!! Creation topography :
432	!! Choice longitudes val1 and val2
433	write(,)'Choice range longitude (in deg.), and topography'
434	ierr=1
435	do while (ierr.ne.0)
436	read(,,iostat=ierr) val1,val2
437	write(,)'Values are:',val1,val2
438	if (ierr.eq.0) then ! got a value
439	! do a sanity check
440	if((val1.lt.-180.).or.(val1.gt.180.)) then
441	write(,)'Longitude should be between -180 and 180 deg!'
442	ierr=1
443	else if((val2.lt.-180.).or.(val2.gt.180.)) then
444	write(,)'Longitude should be between -180 and 180 deg!'
445	ierr=1
446	else if(val1.ge.val2) then
447	write(,)'Longitude 1 should be lower than longitude 2 !'
448	ierr=1
449	else ! find corresponding iref (nearest longitude)
450	! note: rlonv(:) contains increasing values of longitude
451	! starting from -PI to PI
452	do i=1,iip1
453	! get val1
454	write(,) 'val1:',rlonv(i)180./pi,rlonv(i+1)180./pi
455	if ((rlonv(i)*180./pi.le.val1).and.
456	& (rlonv(i+1)*180./pi.ge.val1)) then
457	! find which grid point is nearest to val:
458	if (abs(rlonv(i)*180./pi-val1).le.
459	& abs((rlonv(i+1)*180./pi-val1))) then
460	iref1=i
461	else
462	iref1=i+1
463	endif
464	write(,)'Will use nearest grid longitude 1:',
465	& rlonv(iref1)*180./pi
466	endif
467	enddo ! of do j=1,jjp1
468	do i=1,iip1
469	! get val2
470	if ((rlonv(i)*180./pi.le.val2).and.
471	& (rlonv(i+1)*180./pi.ge.val2)) then
472	! find which grid point is nearest to val:
473	if (abs(rlonv(i)*180./pi-val2).le.
474	& abs((rlonv(i+1)*180./pi-val2))) then
475	iref2=i
476	else
477	iref2=i+1
478	endif
479	write(,)'Will use nearest grid longitude 2:',
480	& rlonv(iref2)*180./pi
481	endif
482	enddo ! of do j=1,jjp1
483	endif ! of if((val.lt.0.).or.(val.gt.90))
484	endif ! of if((val.lt.0.).or.(val.gt.90))
485	enddo ! of do while
486
487	!! Choice latitudes val1 and val2
488	write(,)'Choice range latitudes (in deg.), and topography'
489	ierr=1
490	do while (ierr.ne.0)
491	read(,,iostat=ierr) val1,val2
492	write(,)'Values are:',val1,val2
493	if (ierr.eq.0) then ! got a value
494	! do a sanity check
495	if((val1.lt.-90.).or.(val1.gt.90.)) then
496	write(,)'Latitude should be between -90 and 90 deg!'
497	ierr=1
498	else if((val2.lt.-90.).or.(val2.gt.90.)) then
499	write(,)'Latitude should be between -90 and 90 deg!'
500	ierr=1
501	else if(val1.ge.val2) then
502	write(,)'Latitude 1 should be higuer than latitude 2 !'
503	ierr=1
504	else ! find corresponding jref (nearest latitude)
505	! note: rlatu(:) contains decreasing values of latitude
506	! starting from PI/2 to -PI/2
507	do j=1,jjp1
508	! get val1
509	write(,) 'val1:',rlatu(j)180./pi,rlatu(j+1)180./pi
510	if ((rlatu(j)*180./pi.ge.val1).and.
511	& (rlatu(j+1)*180./pi.le.val1)) then
512	! find which grid point is nearest to val:
513	if (abs(rlatu(j)*180./pi-val1).le.
514	& abs((rlatu(j+1)*180./pi-val1))) then
515	jref1=j
516	else
517	jref1=j+1
518	endif
519	write(,)'Will use nearest grid latitude:',
520	& rlatu(jref1)*180./pi
521	endif
522	enddo ! of do j=1,jjp1
523	do j=1,jjp1
524	! get val2
525	write(,) 'val2:',rlatu(j)180./pi,rlatu(j+1)180./pi
526	if ((rlatu(j)*180./pi.ge.val2).and.
527	& (rlatu(j+1)*180./pi.le.val2)) then
528	! find which grid point is nearest to val:
529	if (abs(rlatu(j)*180./pi-val2).le.
530	& abs((rlatu(j+1)*180./pi-val2))) then
531	jref2=j
532	else
533	jref2=j+1
534	endif
535	write(,)'Will use nearest grid latitude:',
536	& rlatu(jref2)*180./pi
537	endif
538	enddo ! of do j=1,jjp1
539	endif ! of if((val.lt.0.).or.(val.gt.90))
540	endif !
541	enddo ! of do while
542
543	write(,) 'Choice of topography (m) (eg: 1000 or -1000) ?'
544	301 read(,,iostat=ierr) val
545	if(ierr.ne.0) goto 301
546	write(,) 'gravity g is : ',g
547	geop=g*val
548	c phis(:,:)=0
549	phis(iref1:iref2,jref2:jref1)=geop
550	write(,) 'phis=',phis
551	write(,) 'iip1=',iip1,'jjp1=',jjp1
552	write(,) 'iref1=',iref1,'iref2=',iref2
553	write(,) 'jref1=',jref1,'jref2=',jref2
554
555	write(,) 'Do you want another topography choice ?'
556	read(*,fmt='(a)') yes
557	if (yes.eq.'y') goto 302
558	endif ! 'yes flat terrain'
559
560	zmeaS(:,:)=0
561	zstdS(:,:)=0
562	zsigS(:,:)=0
563	zgamS(:,:)=0
564	ztheS(:,:)=0
565
566	yes=' '
567	do while ((yes.ne.'y').and.(yes.ne.'n'))
568	write(,)'Do you want to change soil albedo and IT (y/n)?'
569	read(*,fmt='(a)') yes
570	end do
571	if (yes.eq.'y') then
572	write(,) "Enter value of albedo of the bare ground:"
573	read(,) alb(1,1)
574	alb(:,:)=alb(1,1)
575
576	write(,) "Enter value of thermal inertia of soil:"
577	read(,) surfith(1,1)
578	surfith(:,:)=surfith(1,1)
579
580	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,surfith,surfithfi)
581	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
582	endif !yes for alb and IT changes
583	else
584	yes=' '
585	do while ((yes.ne.'y').and.(yes.ne.'n'))
586	write(,)'Do you want to change soil albedo and IT (y/n)?'
587	read(*,fmt='(a)') yes
588	end do
589	if (yes.eq.'y') then
590	write(,) "Enter value of albedo of the bare ground:"
591	read(,) alb(1,1)
592	alb(:,:)=alb(1,1)
593
594	write(,) "Enter value of thermal inertia of soil:"
595	read(,) surfith(1,1)
596	surfith(:,:)=surfith(1,1)
597
598	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,surfith,surfithfi)
599	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
600	endif !yes for alb and IT changes
601
602	endif ! of if (choice)
603
604	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
605	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zmeaS,zmea)
606	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zstdS,zstd)
607	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zsigS,zsig)
608	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zgamS,zgam)
609	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ztheS,zthe)
610
611	endif ! of if (choix_1.ne.1)
612
613
614	c=======================================================================
615	c Lecture des fichiers (start ou start_archive)
616	c=======================================================================
617
618	if (choix_1.eq.0) then
619
620	write(,) 'Reading file START_ARCHIVE'
621	CALL lect_start_archive(date,tsurf,tsoil,emis,q2,
622	. t,ucov,vcov,ps,teta,phisold_newgrid,q,qsurf,
623	& surfith,nid)
624	write(,) "OK, read start_archive file"
625	! copy soil thermal inertia
626	ithfi(:,:)=inertiedat(:,:)
627
628	ierr= NF_CLOSE(nid)
629
630	else if (choix_1.eq.1) then ! c'est l'appel a tabfi de phyeta0 qui
631	! permet de changer les valeurs du
632	! tab_cntrl Lmodif=1
633	tab0=0
634	Lmodif=1 ! Lmodif set to 1 to allow modifications in phyeta0
635	write(,) 'Reading file START'
636	fichnom = 'start.nc'
637	CALL dynetat0(fichnom,nqmx,vcov,ucov,teta,q,masse,
638	. ps,phis,time)
639
640	write(,) 'Reading file STARTFI'
641	fichnom = 'startfi.nc'
642	CALL phyetat0 (fichnom,tab0,Lmodif,nsoilmx,nqmx,
643	. day_ini,time,
644	. tsurf,tsoil,emis,q2,qsurf)
645
646	! copy albedo and soil thermal inertia
647	do i=1,ngridmx
648	albfi(i) = albedodat(i)
649	do j=1,nsoilmx
650	ithfi(i,j) = inertiedat(i,j)
651	enddo
652	! build a surfithfi(:) using 1st layer of ithfi(:), which might
653	! be neede later on if reinitializing soil thermal inertia
654	surfithfi(i)=ithfi(i,1)
655	enddo
656
657	else
658	CALL exit(1)
659	endif
660
661	dtvr = daysec/FLOAT(day_step)
662	dtphys = dtvr * FLOAT(iphysiq)
663
664	c=======================================================================
665	c
666	c=======================================================================
667
668	do ! infinite loop on list of changes
669
670	write(,)
671	write(,)
672	write(,) 'List of possible changes :'
673	write(,) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
674	write(,)
675	write(,) 'flat : no topography ("aquaplanet")'
676	write(,) 'bilball : uniform albedo and thermal inertia'
677	write(,) 'coldspole : cold subsurface and high albedo at S.pole'
678	write(,) 'qname : change tracer name'
679	write(,) 'q=0 : ALL tracer =zero'
680	write(,) 'q=x : give a specific uniform value to one tracer'
681	write(,) 'ini_q : tracers initialisation for chemistry, water an
682	$d ice '
683	write(,) 'ini_q-H2O : tracers initialisation for chemistry and
684	$ice '
685	write(,) 'ini_q-iceH2O : tracers initialisation for chemistry on
686	$ly '
687	write(,) 'noglacier : Remove tropical H2O ice if \|lat\|<45'
688	write(,) 'nitrogencapn : N2 ice on permanent N polar cap '
689	write(,) 'nitrogencaps : N2 ice on permanent S polar cap '
690	write(,) 'oborealis : H2O ice across Vastitas Borealis'
691	write(,) 'wetstart : start with a wet atmosphere'
692	write(,) 'isotherm : Isothermal Temperatures, wind set to zero'
693	write(,) 'initsoil : initialize soil Temperatures'
694	write(,) 'coldstart : Start X K above the CO2 frost point and
695	$set wind to zero (assumes 100% CO2)'
696	write(,) 'co2ice=0 : remove CO2 polar cap'
697	write(,) 'ptot : change total pressure'
698	write(,) 'emis : change surface emissivity'
699	write(,) 'therm_ini_s: Set soil thermal inertia to reference sur
700	&face values'
701	write(,) 'inert3d: give uniform value of thermal inertial'
702	write(,) 'subsoilice_n: Put deep underground ice layer in northe
703	&rn hemisphere'
704	write(,) 'subsoilice_s: Put deep underground ice layer in southe
705	&rn hemisphere'
706	write(,) 'mons_ice: Put underground ice layer according to MONS-
707	&derived data'
708	write(,) 'plutomap: initialize surface like pluto from ...'
709	write(,) 'subsoil_n2: Put deep underground n2 ice layer where
710	& there is already n2 ice'
711	write(,) 'initsurf: initial surface temperature'
712
713	write(,)
714	write(,) 'Change to perform ?'
715	write(,) ' (enter keyword or return to end)'
716	write(,)
717
718	read(*,fmt='(a20)') modif
719	if (modif(1:1) .eq. ' ') exit ! exit loop on changes
720
721	write(,)
722	write(,) trim(modif) , ' : '
723
724	c 'flat : no topography ("aquaplanet")'
725	c -------------------------------------
726	if (modif(1:len_trim(modif)) .eq. 'flat') then
727	! set topo to zero
728
729	CALL initial0(ip1jmp1,z_reel)
730	CALL multscal(ip1jmp1,z_reel,g,phis)
731	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
732	write(,) 'topography set to zero.'
733
734	c Choice for surface pressure
735	yes=' '
736	do while ((yes.ne.'y').and.(yes.ne.'n'))
737	write(,) 'Do you wish to choose homogeneous surface',
738	& 'pressure (y) or let newstart interpolate ',
739	& ' the previous field (n)?'
740	read(*,fmt='(a)') yes
741	end do
742	if (yes.eq.'y') then
743	write(,) 'New value for ps (Pa) ?'
744	201 read(,,iostat=ierr) patm
745	if(ierr.ne.0) goto 201
746	write(,)
747	write(,) ' new ps everywhere (Pa) = ', patm
748	write(,)
749	do j=1,jjp1
750	do i=1,iip1
751	ps(i,j)=patm
752	enddo
753	enddo
754	end if
755
756	c bilball : albedo, inertie thermique uniforme
757	c --------------------------------------------
758	else if (modif(1:len_trim(modif)) .eq. 'bilball') then
759	write(,) 'constante albedo and iner.therm:'
760	write(,) 'New value for albedo (ex: 0.25) ?'
761	101 read(,,iostat=ierr) alb_bb
762	if(ierr.ne.0) goto 101
763	write(,) ' uniform albedo (new value):',alb_bb
764
765	write(,) 'New value for thermal inertia (eg: 247) ?'
766	150 read(,,iostat=ierr) ith_bb
767	if(ierr.ne.0) goto 150
768	write(,) 'uniform thermal inertia (new value):',ith_bb
769	DO ig=1,ngridmx
770	DO j=1,nsoilmx
771	if(qsurf(ig,igcm_ch4_ice).gt.0)then
772	ithfi(ig,j)=ith_bb
773	else
774	ithfi(ig,j)=400
775	endif
776	ENDDO
777	ENDDO
778	DO j=1,jjp1
779	DO i=1,iip1
780	alb(i,j) = alb_bb ! albedo
781	do isoil=1,nsoilmx
782	ith(i,j,isoil) = ith_bb ! thermal inertia
783	enddo
784	END DO
785	END DO
786	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ith,ithfi)
787	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
788	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
789
790	c TB15 : subsoil_n2 : choice of thermal inertia values for N2 only
791	c ----------------------------------------------------------------
792	else if (modif(1:len_trim(modif)) .eq. 'subsoil_n2') then
793
794	write(,) 'New value for subsoil thermal inertia ?'
795	102 read(,,iostat=ierr) ith_bb
796	if(ierr.ne.0) goto 102
797	write(,) 'thermal inertia (new value):',ith_bb
798
799	write(,)'At which depth (in m.) does the ice layer begin?'
800	write(,)'(currently, the deepest soil layer extends down to:'
801	& ,layer(1),' - ',layer(nsoilmx),')'
802	write(,)'write 0 for uniform value for all subsurf levels?'
803	ierr=1
804	do while (ierr.ne.0)
805	read(,,iostat=ierr) val2
806	write(,)'val2:',val2,'ierr=',ierr
807	if (ierr.eq.0) then ! got a value, but do a sanity check
808	if(val2.gt.layer(nsoilmx)) then
809	write(,)'Depth should be less than ',layer(nsoilmx)
810	ierr=1
811	endif
812	if(val2.lt.layer(1)) then
813	if(val2.eq.0) then
814	write(,)'Thermal inertia value set for all subsurface layers'
815	ierr=0
816	else
817	write(,)'Depth should be more than ',layer(1)
818	ierr=1
819	endif
820	endif
821	endif
822	enddo ! of do while
823
824	! find the reference index iref the depth corresponds to
825	if(val2.eq.0) then
826	iref=1
827	write(,)'Level selected is first level: ',layer(iref),' m'
828	else
829	do isoil=1,nsoilmx-1
830	if((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1)))
831	& then
832	iref=isoil
833	write(,)'Level selected : ',layer(isoil),' m'
834	exit
835	endif
836	enddo
837	endif
838
839	DO ig=1,ngridmx
840	DO j=iref,nsoilmx
841	if((qsurf(ig,igcm_ch4_ice).lt.1.).and.
842	& (qsurf(ig,igcm_co_ice).lt.1.))then
843	ithfi(ig,j)=ith_bb
844	endif
845	ENDDO
846	ENDDO
847
848	CALL gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
849
850	c TB15 Choice tsurf
851	c -----------------------------------------------------
852	else if (modif(1:len_trim(modif)) .eq. 'initsurf') then
853
854	write(,) 'New value for initial surface temperature ?'
855	105 read(,,iostat=ierr) tsurf_bb
856	if(ierr.ne.0) goto 105
857	write(,) 'new surface temperature (K):',tsurf_bb
858	DO ig=1,ngridmx
859	tsurf(ig)=tsurf_bb
860	ENDDO
861
862	c coldspole : sous-sol de la calotte sud toujours froid
863	c -----------------------------------------------------
864	else if (modif(1:len_trim(modif)) .eq. 'coldspole') then
865	write(,)'Sorry, not available'
866
867	c ptot : Modification of the total pressure: ice + current atmosphere
868	c -------------------------------------------------------------------
869	else if (modif(1:len_trim(modif)).eq.'ptot') then
870
871	c calcul de la pression totale glace + atm actuelle
872	patm=0.
873	airetot=0.
874	pcap=0.
875	DO j=1,jjp1
876	DO i=1,iim
877	ig=1+(j-2)*iim +i
878	if(j.eq.1) ig=1
879	if(j.eq.jjp1) ig=ngridmx
880	patm = patm + ps(i,j)*aire(i,j)
881	airetot= airetot + aire(i,j)
882	ENDDO
883	ENDDO
884	ptoto = pcap + patm
885
886	print*,'Current total pressure at surface ',
887	& ptoto/airetot
888
889	print*,'new value?'
890	read(,) ptotn
891	ptotn=ptotn*airetot
892	patmn=ptotn-pcap
893	print*,'ptoto,patm,ptotn,patmn'
894	print*,ptoto,patm,ptotn,patmn
895	print*,'Mult. factor for pressure (atm only)', patmn/patm
896	do j=1,jjp1
897	do i=1,iip1
898	ps(i,j)=ps(i,j)*patmn/patm
899	enddo
900	enddo
901
902	c Correction pour la conservation des traceurs
903	yes=' '
904	do while ((yes.ne.'y').and.(yes.ne.'n'))
905	write(,) 'Do you wish to conserve tracer total mass (y)',
906	& ' or tracer mixing ratio (n) ?'
907	read(*,fmt='(a)') yes
908	end do
909
910	if (yes.eq.'y') then
911	write(,) 'OK : conservation of tracer total mass'
912	DO iq =1, nqmx
913	DO l=1,llm
914	DO j=1,jjp1
915	DO i=1,iip1
916	q(i,j,l,iq)=q(i,j,l,iq)*patm/patmn
917	ENDDO
918	ENDDO
919	ENDDO
920	ENDDO
921	else
922	write(,) 'OK : conservation of tracer mixing ratio'
923	end if
924
925	c qname : change tracer name
926	c --------------------------
927	else if (trim(modif).eq.'qname') then
928	yes='y'
929	do while (yes.eq.'y')
930	write(,) 'Which tracer name do you want to change ?'
931	do iq=1,nqmx
932	write(*,'(i3,a3,a20)')iq,' : ',trim(tnom(iq))
933	enddo
934	write(*,'(a35,i3)')
935	& '(enter tracer number; between 1 and ',nqmx
936	write(,)' or any other value to quit this option)'
937	read(,) iq
938	if ((iq.ge.1).and.(iq.le.nqmx)) then
939	write(,)'Change tracer name ',trim(tnom(iq)),' to ?'
940	read(,) txt
941	tnom(iq)=txt
942	write(,)'Do you want to change another tracer name (y/n)?'
943	read(*,'(a)') yes
944	else
945	! inapropiate value of iq; quit this option
946	yes='n'
947	endif ! of if ((iq.ge.1).and.(iq.le.nqmx))
948	enddo ! of do while (yes.ne.'y')
949
950	c q=0 : set tracers to zero
951	c -------------------------
952	else if (modif(1:len_trim(modif)).eq.'q=0') then
953	c mise a 0 des q (traceurs)
954	write(,) 'Tracers set to 0 (1.E-30 in fact)'
955	DO iq =1, nqmx
956	DO l=1,llm
957	DO j=1,jjp1
958	DO i=1,iip1
959	q(i,j,l,iq)=1.e-30
960	ENDDO
961	ENDDO
962	ENDDO
963	ENDDO
964
965	c set surface tracers to zero
966	DO iq =1, nqmx
967	DO ig=1,ngridmx
968	qsurf(ig,iq)=0.
969	ENDDO
970	ENDDO
971
972	c q=x : initialise tracer manually
973	c --------------------------------
974	else if (modif(1:len_trim(modif)).eq.'q=x') then
975	write(,) 'Which tracer do you want to modify ?'
976	do iq=1,nqmx
977	write(,)iq,' : ',trim(tnom(iq)),' : ',q(1,1,1,iq)
978	enddo
979	write(,) '(choose between 1 and ',nqmx,')'
980	read(,) iq
981	write(,)'mixing ratio of tracer ',trim(tnom(iq)),
982	& ' ? (kg/kg)'
983	read(,) val
984	DO l=1,llm
985	DO j=1,jjp1
986	DO i=1,iip1
987	q(i,j,l,iq)=val
988	ENDDO
989	ENDDO
990	ENDDO
991	yes=' '
992	do while ((yes.ne.'y').and.(yes.ne.'n'))
993	write(,)'Do you want to change surface value (y/n)?'
994	read(*,fmt='(a)') yes
995	end do
996	if (yes.eq.'y') then
997
998	write(,) 'SURFACE value of tracer ',trim(tnom(iq)),
999	& ' ? (kg/m2)'
1000	read(,) val
1001
1002	DO ig=1,ngridmx
1003	qsurf(ig,iq)= val
1004	ENDDO
1005	endif
1006	c DO ig=1,ngridmx
1007	c if (ig.lt.ngridmx*11/32+1) then
1008	c qsurf(ig,iq)=300
1009	c else if (ig.gt.ngridmx*26/32+1) then
1010	c qsurf(ig,iq)=val
1011	c else
1012	c qsurf(ig,iq)=0
1013	c endif
1014	c ENDDO
1015	c endif
1016
1017	c ini_q : Initialize tracers for chemistry
1018	c -----------------------------------------------
1019	else if (modif(1:len_trim(modif)) .eq. 'ini_q') then
1020	write(,) 'Sorry, not available'
1021
1022	c ini_q-H2O : as above exept for the water vapour tracer
1023	c ------------------------------------------------------
1024	else if (modif(1:len_trim(modif)) .eq. 'ini_q-H2O') then
1025	write(,) 'Sorry, not available'
1026
1027	c ini_q-iceH2O : as above exept for ice et H2O
1028	c -----------------------------------------------
1029	else if (modif(1:len_trim(modif)) .eq. 'ini_q-iceH2O') then
1030	write(,) 'Sorry, not available'
1031
1032	c wetstart : wet atmosphere with a north to south gradient
1033	c --------------------------------------------------------
1034	else if (modif(1:len_trim(modif)) .eq. 'wetstart') then
1035	write(,) "No water vapour tracer! Can't use this option"
1036
1037	c noglacier : remove tropical water ice (to initialize high res sim)
1038	c --------------------------------------------------
1039	else if (modif(1:len_trim(modif)) .eq. 'noglacier') then
1040	write(,) "No water vapour tracer! Can't use this option"
1041
1042	c watercapn : H20 ice on permanent northern cap
1043	c --------------------------------------------------
1044	else if (modif(1:len_trim(modif)) .eq. 'watercapn') then
1045	write(,) "No water vapour tracer! Can't use this option"
1046
1047	c watercaps : H20 ice on permanent southern cap
1048	c -------------------------------------------------
1049	else if (modif(1:len_trim(modif)) .eq. 'watercaps') then
1050	write(,) "No water vapour tracer! Can't use this option"
1051
1052	c TB14 surface emissivity
1053	c --------------------
1054	else if (modif(1:len_trim(modif)) .eq. 'emis') then
1055	write(,) 'Emissivity actual max/min values:'
1056	& ,minval(emis),maxval(emis)
1057	write(,) 'Values are in set in physiq !'
1058	write(,) 'Value emissivity N2'
1059	read(,) val
1060	write(,) 'Value emissivity ch4 '
1061	read(,) val2
1062	write(,) 'Value emissivity tholins '
1063	read(,) val3
1064	write(,) 'Value emissivity co '
1065	read(,) val4
1066	write(,) 'New emissivity values :'
1067	write(,) 'N2= ',val,' CH4= ',val2
1068	write(,)' tholins= ',val3,'co= ',val4
1069	DO ig=1,ngridmx
1070	if ((qsurf(ig,igcm_n2).eq.0).and.
1071	& (qsurf(ig,igcm_ch4_ice).eq.0).and.
1072	& (qsurf(ig,igcm_co_ice).eq.0)) then
1073	emis(ig)=val3
1074	else if ((qsurf(ig,igcm_n2).eq.0).and.
1075	& (qsurf(ig,igcm_ch4_ice).eq.0)) then
1076	emis(ig)=val4
1077	else if ((qsurf(ig,igcm_n2).eq.0).and.
1078	& (qsurf(ig,igcm_co_ice).eq.0)) then
1079	emis(ig)=val2
1080	else
1081	emis(ig)=val
1082	endif
1083	ENDDO
1084
1085	c isothermal temperatures and no winds
1086	c ------------------------------------------------
1087	else if (modif(1:len_trim(modif)) .eq. 'isotherm') then
1088	write(,)'Sorry, not available'
1089
1090	c initsoil: subsurface temperature
1091	c ---------------------------------
1092	else if (modif(1:len_trim(modif)) .eq. 'initsoil') then
1093
1094	write(,)'Isothermal temperature of the subsurface'
1095	write(,) 'Value of this temperature ? :'
1096	303 read(,,iostat=ierr) Tiso
1097	if(ierr.ne.0) goto 303
1098
1099	do l=1,nsoilmx
1100	do ig=1, ngridmx
1101	tsoil(ig,l) = Tiso
1102	end do
1103	end do
1104
1105	c -----------------------------------------------------------------------
1106	c coldstart : T set 1K above CO2 frost point and no winds
1107	c ------------------------------------------------
1108	else if (modif(1:len_trim(modif)) .eq. 'coldstart') then
1109	write(,)'Sorry, not available : No CO2 tracer'
1110
1111	c n2ice=0 : remove n2 polar ice caps'
1112	c ------------------------------------------------
1113	else if (modif(1:len_trim(modif)) .eq. 'N2ice=0') then
1114	! check that there is indeed a N2_ice tracer ...
1115	if (igcm_n2.ne.0) then
1116	do ig=1,ngridmx
1117	!co2ice(ig)=0
1118	qsurf(ig,igcm_n2)=0
1119	emis(ig)=emis(ngridmx/2)
1120	end do
1121	else
1122	write(,) "Can't remove N2 ice!! (no N2 tracer)"
1123	endif
1124
1125	! therm_ini_s: (re)-set soil thermal inertia to reference surface values
1126	! ----------------------------------------------------------------------
1127	else if (modif(1:len_trim(modif)).eq.'therm_ini_s') then
1128	! write(,)"surfithfi(1):",surfithfi(1)
1129	do isoil=1,nsoilmx
1130	inertiedat(1:ngridmx,isoil)=surfithfi(1:ngridmx)
1131	enddo
1132	write(,)'OK: Soil thermal inertia has been reset to referenc
1133	&e surface values'
1134	write(,)"inertiedat(1,1):",inertiedat(1,1)
1135	ithfi(:,:)=inertiedat(:,:)
1136	! recast ithfi() onto ith()
1137	call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
1138	! Check:
1139	! do i=1,iip1
1140	! do j=1,jjp1
1141	! do isoil=1,nsoilmx
1142	! write(77,*) i,j,isoil," ",ith(i,j,isoil)
1143	! enddo
1144	! enddo
1145	! enddo
1146
1147	! TB14 inert3d: set soil thermal inertia to specific uniform value
1148	! ----------------------------------------------------------------------
1149	else if (modif(1:len_trim(modif)).eq.'inert3d') then
1150	write(,) 'Actual value of surf Thermal inertia at ig=1: ',
1151	& inertiedat(1,1), ' SI'
1152	write(,) 'Value of Thermal inertia (SI) ? '
1153	read(,) val
1154	do isoil=1,nsoilmx
1155	do ig=1,ngridmx
1156	inertiedat(ig,isoil)=val
1157	enddo
1158	enddo
1159	write(,)'OK: Soil TI set to ',val,' SI everywhere'
1160	ithfi(:,:)=inertiedat(:,:)
1161	! recast ithfi() onto ith()
1162	call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
1163
1164	! subsoilice_n: Put deep ice layer in northern hemisphere soil
1165	! ------------------------------------------------------------
1166	else if (modif(1:len_trim(modif)).eq.'subsoilice_n') then
1167
1168	write(,)'From which latitude (in deg.), up to the north pole,
1169	&should we put subterranean ice?'
1170	ierr=1
1171	do while (ierr.ne.0)
1172	read(,,iostat=ierr) val
1173	if (ierr.eq.0) then ! got a value
1174	! do a sanity check
1175	if((val.lt.0.).or.(val.gt.90)) then
1176	write(,)'Latitude should be between 0 and 90 deg. !!!'
1177	ierr=1
1178	else ! find corresponding jref (nearest latitude)
1179	! note: rlatu(:) contains decreasing values of latitude
1180	! starting from PI/2 to -PI/2
1181	do j=1,jjp1
1182	if ((rlatu(j)*180./pi.ge.val).and.
1183	& (rlatu(j+1)*180./pi.le.val)) then
1184	! find which grid point is nearest to val:
1185	if (abs(rlatu(j)*180./pi-val).le.
1186	& abs((rlatu(j+1)*180./pi-val))) then
1187	jref=j
1188	else
1189	jref=j+1
1190	endif
1191
1192	write(,)'Will use nearest grid latitude which is:',
1193	& rlatu(jref)*180./pi
1194	endif
1195	enddo ! of do j=1,jjp1
1196	endif ! of if((val.lt.0.).or.(val.gt.90))
1197	endif !of if (ierr.eq.0)
1198	enddo ! of do while
1199
1200	! Build layers() (as in soil_settings.F)
1201	val2=sqrt(mlayer(0)*mlayer(1))
1202	val3=mlayer(1)/mlayer(0)
1203	do isoil=1,nsoilmx
1204	layer(isoil)=val2(val3*(isoil-1))
1205	enddo
1206
1207	write(,)'At which depth (in m.) does the ice layer begin?'
1208	write(,)'(currently, the deepest soil layer extends down to:'
1209	& ,layer(nsoilmx),')'
1210	ierr=1
1211	do while (ierr.ne.0)
1212	read(,,iostat=ierr) val2
1213	! write(,)'val2:',val2,'ierr=',ierr
1214	if (ierr.eq.0) then ! got a value, but do a sanity check
1215	if(val2.gt.layer(nsoilmx)) then
1216	write(,)'Depth should be less than ',layer(nsoilmx)
1217	ierr=1
1218	endif
1219	if(val2.lt.layer(1)) then
1220	write(,)'Depth should be more than ',layer(1)
1221	ierr=1
1222	endif
1223	endif
1224	enddo ! of do while
1225
1226	! find the reference index iref the depth corresponds to
1227	! if (val2.lt.layer(1)) then
1228	! iref=1
1229	! else
1230	do isoil=1,nsoilmx-1
1231	if((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1)))
1232	& then
1233	iref=isoil
1234	exit
1235	endif
1236	enddo
1237	! endif
1238
1239	! write(,)'iref:',iref,' jref:',jref
1240	! write(,)'layer',layer
1241	! write(,)'mlayer',mlayer
1242
1243	! thermal inertia of the ice:
1244	ierr=1
1245	do while (ierr.ne.0)
1246	write(,)'What is the value of subterranean ice thermal inert
1247	&ia? (e.g.: 2000)'
1248	read(,,iostat=ierr)iceith
1249	enddo ! of do while
1250
1251	! recast ithfi() onto ith()
1252	call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
1253
1254	do j=1,jref
1255	! write(,)'j:',j,'rlatu(j)180./pi:',rlatu(j)180./pi
1256	do i=1,iip1 ! loop on longitudes
1257	! Build "equivalent" thermal inertia for the mixed layer
1258	ith(i,j,iref+1)=sqrt((layer(iref+1)-layer(iref))/
1259	& (((val2-layer(iref))/(ith(i,j,iref)**2))+
1260	& ((layer(iref+1)-val2)/(iceith)**2)))
1261	! Set thermal inertia of lower layers
1262	do isoil=iref+2,nsoilmx
1263	ith(i,j,isoil)=iceith ! ice
1264	enddo
1265	enddo ! of do i=1,iip1
1266	enddo ! of do j=1,jjp1
1267
1268
1269	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
1270
1271	! do i=1,nsoilmx
1272	! write(,)'i:',i,'ithfi(1,i):',ithfi(1,i)
1273	! enddo
1274
1275
1276	! subsoilice_s: Put deep ice layer in southern hemisphere soil
1277	! ------------------------------------------------------------
1278	else if (modif(1:len_trim(modif)).eq.'subsoilice_s') then
1279
1280	write(,)'From which latitude (in deg.), down to the south pol
1281	&e, should we put subterranean ice?'
1282	ierr=1
1283	do while (ierr.ne.0)
1284	read(,,iostat=ierr) val
1285	if (ierr.eq.0) then ! got a value
1286	! do a sanity check
1287	if((val.gt.0.).or.(val.lt.-90)) then
1288	write(,)'Latitude should be between 0 and -90 deg. !!!'
1289	ierr=1
1290	else ! find corresponding jref (nearest latitude)
1291	! note: rlatu(:) contains decreasing values of latitude
1292	! starting from PI/2 to -PI/2
1293	do j=1,jjp1
1294	if ((rlatu(j)*180./pi.ge.val).and.
1295	& (rlatu(j+1)*180./pi.le.val)) then
1296	! find which grid point is nearest to val:
1297	if (abs(rlatu(j)*180./pi-val).le.
1298	& abs((rlatu(j+1)*180./pi-val))) then
1299	jref=j
1300	else
1301	jref=j+1
1302	endif
1303
1304	write(,)'Will use nearest grid latitude which is:',
1305	& rlatu(jref)*180./pi
1306	endif
1307	enddo ! of do j=1,jjp1
1308	endif ! of if((val.lt.0.).or.(val.gt.90))
1309	endif !of if (ierr.eq.0)
1310	enddo ! of do while
1311
1312	! Build layers() (as in soil_settings.F)
1313	val2=sqrt(mlayer(0)*mlayer(1))
1314	val3=mlayer(1)/mlayer(0)
1315	do isoil=1,nsoilmx
1316	layer(isoil)=val2(val3*(isoil-1))
1317	enddo
1318
1319	write(,)'At which depth (in m.) does the ice layer begin?'
1320	write(,)'(currently, the deepest soil layer extends down to:'
1321	& ,layer(nsoilmx),')'
1322	ierr=1
1323	do while (ierr.ne.0)
1324	read(,,iostat=ierr) val2
1325	! write(,)'val2:',val2,'ierr=',ierr
1326	if (ierr.eq.0) then ! got a value, but do a sanity check
1327	if(val2.gt.layer(nsoilmx)) then
1328	write(,)'Depth should be less than ',layer(nsoilmx)
1329	ierr=1
1330	endif
1331	if(val2.lt.layer(1)) then
1332	write(,)'Depth should be more than ',layer(1)
1333	ierr=1
1334	endif
1335	endif
1336	enddo ! of do while
1337
1338	! find the reference index iref the depth corresponds to
1339	do isoil=1,nsoilmx-1
1340	if((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1)))
1341	& then
1342	iref=isoil
1343	exit
1344	endif
1345	enddo
1346
1347	! thermal inertia of the ice:
1348	ierr=1
1349	do while (ierr.ne.0)
1350	write(,)'What is the value of subterranean ice thermal inert
1351	&ia? (e.g.: 2000)'
1352	read(,,iostat=ierr)iceith
1353	enddo ! of do while
1354
1355	! recast ithfi() onto ith()
1356	call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
1357
1358	do j=jref,jjp1
1359	! write(,)'j:',j,'rlatu(j)180./pi:',rlatu(j)180./pi
1360	do i=1,iip1 ! loop on longitudes
1361	! Build "equivalent" thermal inertia for the mixed layer
1362	ith(i,j,iref+1)=sqrt((layer(iref+1)-layer(iref))/
1363	& (((val2-layer(iref))/(ith(i,j,iref)**2))+
1364	& ((layer(iref+1)-val2)/(iceith)**2)))
1365	! Set thermal inertia of lower layers
1366	do isoil=iref+2,nsoilmx
1367	ith(i,j,isoil)=iceith ! ice
1368	enddo
1369	enddo ! of do i=1,iip1
1370	enddo ! of do j=jref,jjp1
1371
1372
1373	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
1374
1375	c 'mons_ice' : use MONS data to build subsurface ice table
1376	c --------------------------------------------------------
1377	else if (modif(1:len_trim(modif)).eq.'mons_ice') then
1378	write(,) "Sorry, not available"
1379
1380	c --------------------------------------------------------
1381	c 'plutomap' : initialize pluto
1382	c --------------------------------------------------------
1383	else if (modif(1:len_trim(modif)).eq.'plutomap') then
1384
1385	write(,) 'Pluto_map.dat has to be in your simulation file !!'
1386	write(,) 'Initial value of CH4 ice (kg)'
1387	read(,,iostat=ierr) val
1388	write(,) 'Initial value of N2 ice (kg)'
1389	read(,,iostat=ierr) val2
1390	write(,) 'Initial value of CO ice (kg)'
1391	read(,,iostat=ierr) val3
1392
1393	c Reading map codes (from Lellouch et al. 2011)
1394	c The map must be placed in your simulation file
1395	open(27,file='pluto_map.dat')
1396	do j=1,jm_plu+1
1397	read(27,*) (map_pluto_dat(i,j),i=1,im_plu+1)
1398	do i=1,im_plu+1
1399	N2_pluto_dat(i,j) =0.
1400	CH4_pluto_dat(i,j) =0.
1401	CO_pluto_dat(i,j) =0.
1402	if(map_pluto_dat(i,j).eq.1) N2_pluto_dat(i,j)=val2
1403	if(map_pluto_dat(i,j).eq.2) CH4_pluto_dat(i,j)=val
1404	c if(map_pluto_dat(i,j).eq.3) -> nothing : Tholins
1405	if(map_pluto_dat(i,j).eq.4) CO_pluto_dat(i,j)=val3
1406	c if(map_pluto_dat(i,j).eq.1) alb_pluto_dat(i,j)=0.657
1407	c if(map_pluto_dat(i,j).eq.1) ith_pluto_dat(i,j)=100
1408	end do
1409	end do
1410	close (27)
1411
1412	c Option to remove N2 ice in the south as in Young/Olkin model
1413
1414	yes=' '
1415	do while ((yes.ne.'y').and.(yes.ne.'n'))
1416	write(,)'No nitrogen in the south? as in Olkin et al. 2014'
1417	read(*,fmt='(a)') yes
1418	end do
1419	if (yes.eq.'y') then
1420	write(,) 'Lower latitude of nitrogen (deg: -90 to 90)'
1421	read(,) val
1422	itau = nint(jm_plu*(90.-val)/180.)+1
1423	write(,) "grid lower latitude="
1424	write(,)90.-180.*(itau-1)/(jm_plu)
1425	do j=itau, jm_plu+1
1426	do i=1,im_plu+1
1427	if(map_pluto_dat(i,j).eq.1) N2_pluto_dat(i,j)=0.
1428	end do
1429	end do
1430	end if
1431
1432	c Reindexation to shift the longitude grid like a LMD GCM grid...
1433	do j=1,jm_plu+1
1434	N2_pluto_rein(1,j)=(N2_pluto_dat(1,j)+N2_pluto_dat(im_plu,j))/2
1435	CH4_pluto_rein(1,j)=(CH4_pluto_dat(1,j)+CH4_pluto_dat(im_plu,j))/2
1436	CO_pluto_rein(1,j)=(CO_pluto_dat(1,j)+CO_pluto_dat(im_plu,j))/2
1437	c alb_pluto_rein(1,j)=(alb_pluto_dat(1,j)+alb_pluto_dat(im_plu,j))/2
1438
1439	c ith_pluto_rein(1,j)=(ith_pluto_dat(1,j)+ith_pluto_dat(im_plu,j))/2
1440	do i=2,im_plu
1441	N2_pluto_rein(i,j)=(N2_pluto_dat(i-1,j)+N2_pluto_dat(i,j))/2
1442	CH4_pluto_rein(i,j)=(CH4_pluto_dat(i-1,j)+CH4_pluto_dat(i,j))/2
1443	CO_pluto_rein(i,j)=(CO_pluto_dat(i-1,j)+CO_pluto_dat(i,j))/2
1444	c alb_pluto_rein(i,j)=(alb_pluto_dat(i-1,j)+alb_pluto_dat(i,j))/2
1445	c ith_pluto_rein(i,j)=(ith_pluto_dat(i-1,j)+ith_pluto_dat(i,j))/2
1446	end do
1447	N2_pluto_rein(im_plu+1,j) = N2_pluto_rein(1,j)
1448	CH4_pluto_rein(im_plu+1,j) = CH4_pluto_rein(1,j)
1449	CO_pluto_rein(im_plu+1,j) = CO_pluto_rein(1,j)
1450	c alb_pluto_rein(im_plu+1,j) = alb_pluto_rein(1,j)
1451	c ith_pluto_rein(im_plu+1,j)= ith_pluto_dat(1,j)
1452	end do
1453
1454	c latitude and longitude in REindexed pluto map :
1455	latv_plu(1)=90. *pi/180.
1456	do j=2,jm_plu
1457	latv_plu(j)= (90-(j-1 -0.5)180./(jm_plu-1))pi/180.
1458	end do
1459	latv_plu(jm_plu+1)= -90. *pi/180.
1460
1461	do i=1,im_plu+1
1462	lonu_plu(i)=(-180+ (i-1)360./(im_plu) ) pi/180.
1463	enddo
1464	c write(,) 'TB15 im_plu = ',im_plu
1465	c write(,) 'TB15 jm_plu = ',jm_plu
1466	c write(,) 'TB15 imn = ',iim
1467	c write(,) 'TB15 jmn = ',jjm
1468	c write(,) 'TB15 lonu_plu = ',lonu_plu*180./pi
1469	c write(,) 'TB15 latv_plu = ',latv_plu*180./pi
1470	c write(,) 'TB15 rlonu = ',rlonu*180./pi
1471	c write(,) 'TB15 rlatv = ',rlatv*180./pi
1472	c interpolation in LMD GCM grid:
1473	c write(,) 'TB15 N2_pluto_rein = ',N2_pluto_rein
1474	call interp_horiz(N2_pluto_rein,N2_pluto_gcm,
1475	& im_plu,jm_plu,iim,jjm,1,lonu_plu,latv_plu,rlonu,rlatv)
1476	call interp_horiz(CH4_pluto_rein,CH4_pluto_gcm,
1477	& im_plu,jm_plu,iim,jjm,1,lonu_plu,latv_plu,rlonu,rlatv)
1478	call interp_horiz(CO_pluto_rein,CO_pluto_gcm,
1479	& im_plu,jm_plu,iim,jjm,1,lonu_plu,latv_plu,rlonu,rlatv)
1480	c call interp_horiz(alb_pluto_rein,alb_pluto_gcm,
1481	c & im_plu,jm_plu,iim,jjm,1,lonu_plu,latv_plu,rlonu,rlatv)
1482	c write(,) 'TB15 N2_pluto_gcm = ',N2_pluto_gcm
1483	c do i=1,iip1
1484	c do j=1,jjp1
1485	c N2_pluto_gcm(i,j) = 1000.
1486	c end do
1487	c end do
1488	c write(,) 'TB15 N2_pluto_gcm = ',N2_pluto_gcm
1489
1490	c Option to remove N2 ice in the North
1491	c yes=' '
1492	c do while ((yes.ne.'y').and.(yes.ne.'n'))
1493	c write(,)'No N2 ice at North pole ?'
1494	c read(*,fmt='(a)') yes
1495	c end do
1496	c if (yes.eq.'y') then
1497	c do i=1,iip1
1498	c N2_pluto_gcm(i,1) = 0.
1499	c end do
1500	c end if
1501
1502	cc call interp_horiz(ith_pluto_rein,ith_pluto_gcm,
1503	cc & im_plu,jm_plu,iim,jjm,1,lonu_plu,latv_plu,rlonu,rlatv)
1504
1505	qsurf(:,:) =0.
1506	c CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb_pluto_gcm,albfi)
1507	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,
1508	& N2_pluto_gcm,qsurf(1,igcm_n2))
1509	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,
1510	& CH4_pluto_gcm,qsurf(1,igcm_ch4_ice))
1511	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,
1512	& CO_pluto_gcm,qsurf(1,igcm_co_ice))
1513
1514	c ----------------------------------------------------------------
1515	c ----------------------------------------------------------------
1516
1517	else
1518	write(,) ' Unknown (misspelled?) option!!!'
1519	end if ! of if (modif(1:len_trim(modif)) .eq. '...') elseif ...
1520
1521	enddo ! of do ! infinite loop on liste of changes
1522
1523	999 continue
1524
1525
1526	c=======================================================================
1527	c Correct pressure on the new grid (menu 0)
1528	c=======================================================================
1529
1530	if ((choix_1.eq.0).and.(.not.flagps0)) then
1531
1532	r = 8.314511E+0 *1000.E+0/mugaz
1533	do j=1,jjp1
1534	do i=1,iip1
1535	ps(i,j) = ps(i,j) *
1536	. exp((phisold_newgrid(i,j)-phis(i,j)) /
1537	. (t(i,j,1) * r))
1538	end do
1539	end do
1540
1541	c periodicity of surface ps in longitude
1542	do j=1,jjp1
1543	ps(1,j) = ps(iip1,j)
1544	end do
1545	end if
1546
1547	c=======================================================================
1548	c=======================================================================
1549
1550	c=======================================================================
1551	c Initialisation de la physique / ecriture de newstartfi :
1552	c=======================================================================
1553
1554
1555	CALL inifilr
1556	CALL pression(ip1jmp1, ap, bp, ps, p3d)
1557
1558	c-----------------------------------------------------------------------
1559	c Initialisation pks:
1560	c-----------------------------------------------------------------------
1561
1562	CALL exner_hyb(ip1jmp1, ps, p3d, beta, pks, pk, pkf)
1563
1564	! Calcul de la temperature potentielle teta
1565
1566	if (flagtset) then
1567	DO l=1,llm
1568	DO j=1,jjp1
1569	DO i=1,iim
1570	teta(i,j,l) = Tset(i,j,l) * cpp/pk(i,j,l)
1571	ENDDO
1572	teta (iip1,j,l)= teta (1,j,l)
1573	ENDDO
1574	ENDDO
1575	else if (choix_1.eq.0) then
1576	DO l=1,llm
1577	DO j=1,jjp1
1578	DO i=1,iim
1579	teta(i,j,l) = t(i,j,l) * cpp/pk(i,j,l)
1580	ENDDO
1581	teta (iip1,j,l)= teta (1,j,l)
1582	ENDDO
1583	ENDDO
1584	endif
1585
1586	C Calcul intermediaire
1587
1588	if (choix_1.eq.0) then
1589	CALL massdair( p3d, masse )
1590
1591	print *,' ALPHAX ',alphax
1592	DO l = 1, llm
1593	DO i = 1, iim
1594	xppn(i) = aire( i, 1 ) * masse( i , 1 , l )
1595	xpps(i) = aire( i,jjp1 ) * masse( i , jjp1 , l )
1596	ENDDO
1597	xpn = SUM(xppn)/apoln
1598	xps = SUM(xpps)/apols
1599	DO i = 1, iip1
1600	masse( i , 1 , l ) = xpn
1601	masse( i , jjp1 , l ) = xps
1602	ENDDO
1603	ENDDO
1604	endif
1605
1606
1607	phis(iip1,:) = phis(1,:)
1608
1609	CALL inidissip ( lstardis, nitergdiv, nitergrot, niterh,
1610	* tetagdiv, tetagrot , tetatemp )
1611	itau=0
1612
1613	if (choix_1.eq.0) then
1614	day_ini=int(date)
1615	endif
1616
1617	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
1618
1619	CALL caldyn0( itau,ucov,vcov,teta,ps,masse,pk,phis ,
1620	* phi,w, pbaru,pbarv,day_ini+time )
1621
1622	c CALL caldyn
1623	c $ ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
1624	c $ phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, day_ini )
1625	CALL dynredem0("restart.nc",day_ini,anneeref,phis,nqmx)
1626	CALL dynredem1("restart.nc",0.0,vcov,ucov,teta,q,nqmx,masse,ps)
1627
1628	C Ecriture etat initial physique
1629
1630	call physdem1("restartfi.nc",lonfi,latfi,nsoilmx,nqmx,
1631	. dtphys,float(day_ini),
1632	. time,tsurf,tsoil,emis,q2,qsurf,
1633	. airefi,albfi,ithfi,zmea,zstd,zsig,zgam,zthe)
1634
1635	write(,) '***************************************'
1636	write(,) 'FINAL VALUES IN RESTART: '
1637	write(,) 'check consistency'
1638	write(,) 'to change values in restart use start2archive'
1639	write(,) '***************************************'
1640	write(,) 'radius= ',rad
1641	write(,) 'omeg= ',omeg
1642	write(,) 'gravity= ',g
1643	write(,) 'cpp= ', cpp
1644	write(,) 'mugaz= ',mugaz
1645	write(,) 'daysec=',daysec
1646	write(,) 'kappa=', 8.314511E+0 1000.E+0/(mugazcpp)
1647
1648	c=======================================================================
1649	c Formats
1650	c=======================================================================
1651
1652	1 FORMAT(//10x,'la valeur de im =',i4,2x,'lue sur le fichier de dema
1653	*rrage est differente de la valeur parametree iim =',i4//)
1654	2 FORMAT(//10x,'la valeur de jm =',i4,2x,'lue sur le fichier de dema
1655	*rrage est differente de la valeur parametree jjm =',i4//)
1656	3 FORMAT(//10x,'la valeur de lllm =',i4,2x,'lue sur le fichier demar
1657	*rage est differente de la valeur parametree llm =',i4//)
1658
1659	end
1660
1661	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1662	c subroutine load_MONS_data(MONS_Hdn,MONS_d21)
1663

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

Download in other formats:

Original Format