Context Navigation

newstart.F @ 2942

Last change on this file since 2942 was 2942, checked in by llange, 20 months ago
Mars PCM Add the possibility to use a different thermal inertia (field 'inertiesoil') than inertiedat in the PCM (for paleoclimate studies). By defaut, if there is not inertiesoil, inertiedat is used. Soil_tifeedback still work with inertiedat Newstart adapted, start2archive will be modified by Romain Vandemeulebrouck. LL
File size: 74.6 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	use ioipsl_getincom, only: getin
18	use mod_phys_lmdz_para, only: is_parallel, is_sequential,
19	& is_mpi_root, is_omp_root,
20	& is_master
21	use infotrac, only: infotrac_init, nqtot, tname
22	use tracer_mod, only: noms, mmol,
23	& igcm_dust_number, igcm_dust_mass,
24	& igcm_ccn_number, igcm_ccn_mass,
25	& igcm_h2o_vap, igcm_h2o_ice, igcm_co2,
26	& igcm_hdo_vap, igcm_hdo_ice,
27	& igcm_n2, igcm_ar, igcm_o2, igcm_co,
28	& igcm_o, igcm_h2
29	use surfdat_h, only: phisfi, z0, zmea, zstd, zsig, zgam, zthe,
30	& albedodat, z0_default, qsurf, tsurf,
31	& emis, hmons, summit, base, watercap,
32	& ini_surfdat_h_slope_var,end_surfdat_h_slope_var
33	use comsoil_h, only: inertiedat, inertiesoil,layer, mlayer,
34	& nsoilmx, tsoil,ini_comsoil_h_slope_var, end_comsoil_h_slope_var
35	use control_mod, only: day_step, iphysiq, anneeref, planet_type
36	use geometry_mod, only: longitude,latitude,cell_area
37	use phyetat0_mod, only: phyetat0
38	use phyredem, only: physdem0, physdem1
39	use iostart, only: open_startphy
40	use dimradmars_mod, only: albedo,
41	& ini_dimradmars_mod_slope_var,end_dimradmars_mod_slope_var
42	use dust_param_mod, only: tauscaling
43	use turb_mod, only: q2, wstar
44	use filtreg_mod, only: inifilr
45	USE mod_const_mpi, ONLY: COMM_LMDZ
46	USE comvert_mod, ONLY: ap,bp,pa,preff
47	USE comconst_mod, ONLY: lllm,daysec,dtphys,dtvr,
48	. cpp,kappa,rad,omeg,g,r,pi
49	USE serre_mod, ONLY: alphax
50	USE temps_mod, ONLY: day_ini,hour_ini
51	USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
52	USE iniphysiq_mod, ONLY: iniphysiq
53	USE exner_hyb_m, ONLY: exner_hyb
54	USE inichim_newstart_mod, ONLY: inichim_newstart
55	use comslope_mod, ONLY: nslope,def_slope,def_slope_mean,
56	& subslope_dist,end_comslope_h,ini_comslope_h
57	USE subslope_mola_mod, ONLY: subslope_mola
58
59	implicit none
60
61	include "dimensions.h"
62	integer, parameter :: ngridmx = (2+(jjm-1)*iim - 1/jjm)
63	include "paramet.h"
64	include "comgeom2.h"
65	include "comdissnew.h"
66	include "clesph0.h"
67	include "netcdf.inc"
68	c=======================================================================
69	c Declarations
70	c=======================================================================
71
72	c Variables dimension du fichier "start_archive"
73	c------------------------------------
74	CHARACTER relief*3
75
76	c et autres:
77	c----------
78
79	c Variables pour les lectures NetCDF des fichiers "start_archive"
80	c--------------------------------------------------
81	INTEGER nid_dyn, nid_fi,nid,nvarid
82	INTEGER tab0
83
84	REAL date
85	REAL p_rad,p_omeg,p_g,p_mugaz,p_daysec
86
87	c Variable histoire
88	c------------------
89	REAL vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm) ! vents covariants
90	REAL phis(iip1,jjp1)
91	REAL,ALLOCATABLE :: q(:,:,:,:) ! champs advectes
92
93	c autre variables dynamique nouvelle grille
94	c------------------------------------------
95	REAL pks(iip1,jjp1)
96	REAL w(iip1,jjp1,llm+1)
97	REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
98	! REAL dv(ip1jm,llm),du(ip1jmp1,llm)
99	! REAL dh(ip1jmp1,llm),dp(ip1jmp1)
100	REAL phi(iip1,jjp1,llm)
101
102	integer klatdat,klongdat
103	PARAMETER (klatdat=180,klongdat=360)
104
105	c Physique sur grille scalaire
106	c----------------------------
107	real zmeaS(iip1,jjp1),zstdS(iip1,jjp1)
108	real zsigS(iip1,jjp1),zgamS(iip1,jjp1),ztheS(iip1,jjp1)
109	real hmonsS(iip1,jjp1)
110	real summitS(iip1,jjp1)
111	real baseS(iip1,jjp1)
112	real zavgS(iip1,jjp1)
113	real z0S(iip1,jjp1)
114
115	c variable physique
116	c------------------
117	REAL tauscadyn(iip1,jjp1) ! dust conversion factor on the dynamics grid
118	real alb(iip1,jjp1),albfi(ngridmx) ! albedos
119	real ith(iip1,jjp1,nsoilmx),ithfi(ngridmx,nsoilmx) ! thermal inertia (3D)
120	real surfith(iip1,jjp1),surfithfi(ngridmx) ! surface thermal inertia (2D)
121	! REAL latfi(ngridmx),lonfi(ngridmx),airefi(ngridmx)
122
123	INTEGER i,j,l,isoil,ig,idum
124	real mugaz ! molar mass of the atmosphere
125
126	integer ierr !, nbetat
127
128	c Variables on the new grid along scalar points
129	c------------------------------------------------------
130	! REAL p(iip1,jjp1)
131	REAL t(iip1,jjp1,llm)
132	real phisold_newgrid(iip1,jjp1)
133	REAL :: teta(iip1, jjp1, llm)
134	REAL :: pk(iip1,jjp1,llm)
135	REAL :: pkf(iip1,jjp1,llm)
136	REAL :: ps(iip1, jjp1)
137	REAL :: masse(iip1,jjp1,llm)
138	REAL :: xpn,xps,xppn(iim),xpps(iim)
139	REAL :: p3d(iip1, jjp1, llm+1)
140	! REAL dteta(ip1jmp1,llm)
141
142	c Variable de l'ancienne grille
143	c------------------------------
144	real time
145	real tab_cntrl(100)
146	real tab_cntrl_bis(100)
147
148	c variables diverses
149	c-------------------
150	real choix_1 ! ==0 : read start_archive file ; ==1: read start files
151	character*80 fichnom
152	integer Lmodif,iq
153	integer flagthermo, flagh2o
154	character modif*20
155	real tsud,albsud,alb_bb,ith_bb,Tiso
156	real ptoto,pcap,patm,airetot,ptotn,patmn
157	! real ssum
158	character*1 yes
159	logical :: flagiso=.false. , flagps0=.false.
160	real val, val2, val3 ! to store temporary variables
161	real :: iceith=2000 ! thermal inertia of subterranean ice
162	real :: iceithN,iceithS ! values of thermal inertias in N & S hemispheres
163	integer iref,jref
164
165	INTEGER :: itau
166	real DoverH !D/H ratio
167
168	INTEGER :: numvanle
169	character(len=50) :: txt ! to store some text
170	integer :: count
171	real :: profile(llm+1) ! to store an atmospheric profile + surface value
172
173	! MONS data:
174	real :: MONS_Hdn(iip1,jjp1) ! Hdn: %WEH=Mass fraction of H2O
175	real :: MONS_d21(iip1,jjp1) ! ice table "depth" (in kg/m2)
176	! coefficient to apply to convert d21 to 'true' depth (m)
177	real :: MONS_coeff
178	real :: MONS_coeffS ! coeff for southern hemisphere
179	real :: MONS_coeffN ! coeff for northern hemisphere
180	! real,parameter :: icedepthmin=1.e-3 ! Ice begins at most at that depth
181	! Reference position for composition
182	real :: latref,lonref,dlatmin,dlonmin
183	! Variable used to change composition
184	real :: Svmr,Smmr,Smmr_old,Smmr_new,n,Sn
185	real :: Mair_old,Mair_new,vmr_old,vmr_new
186	real,allocatable :: coefvmr(:) ! Correction coefficient when changing composition
187	real :: maxq
188	integer :: iloc(1), iqmax, islope
189	! sub-grid cloud fraction
190	real :: totcloudfrac(ngridmx)
191	!Variable to change the number of subslope
192	REAL, ALLOCATABLE :: default_def_slope(:)
193	REAL,ALLOCATABLE :: tsurf_old_slope(:,:) ! Surface temperature (K)
194	REAL,ALLOCATABLE :: emis_old_slope(:,:) ! Thermal IR surface emissivity
195	REAL,ALLOCATABLE :: qsurf_old_slope(:,:,:) ! tracer on surface (e.g. kg.m-2)
196	REAL,ALLOCATABLE :: watercap_old_slope(:,:) ! Surface water ice (kg.m-2)
197	REAL,ALLOCATABLE :: tsoil_old_slope(:,:,:)
198	REAL,ALLOCATABLE :: albedo_old_slope(:,:,:) ! Surface albedo in each solar band
199	integer :: iflat
200	integer :: nslope_old, nslope_new
201
202
203
204	c sortie visu pour les champs dynamiques
205	c---------------------------------------
206	! INTEGER :: visuid
207	! real :: time_step,t_ops,t_wrt
208	! CHARACTER*80 :: visu_file
209
210	cpp = 744.499 ! for Mars, instead of 1004.70885 (Earth)
211	preff = 610. ! for Mars, instead of 101325. (Earth)
212	pa= 20 ! for Mars, instead of 500 (Earth)
213	planet_type="mars"
214
215	! initialize "serial/parallel" related stuff:
216	! (required because we call tabfi() below, before calling iniphysiq)
217	is_sequential=.true.
218	is_parallel=.false.
219	is_mpi_root=.true.
220	is_omp_root=.true.
221	is_master=.true.
222
223	! Load tracer number and names:
224	call infotrac_init
225	! allocate arrays
226	allocate(q(iip1,jjp1,llm,nqtot))
227	allocate(coefvmr(nqtot))
228
229	c=======================================================================
230	c Choice of the start file(s) to use
231	c=======================================================================
232
233	write(,) 'From which kind of files do you want to create new',
234	. 'start and startfi files'
235	write(,) ' 0 - from a file start_archive'
236	write(,) ' 1 - from files start and startfi'
237
238	c-----------------------------------------------------------------------
239	c Open file(s) to modify (start or start_archive)
240	c-----------------------------------------------------------------------
241
242	DO
243	read(,,iostat=ierr) choix_1
244	if ((choix_1 /= 0).OR.(choix_1 /=1)) EXIT
245	ENDDO
246
247	c Open start_archive
248	c ~~~~~~~~~~~~~~~~~~~~~~~~~~
249	if (choix_1.eq.0) then
250
251	write(,) 'Creating start files from:'
252	write(,) './start_archive.nc'
253	write(,)
254	fichnom = 'start_archive.nc'
255	ierr = NF_OPEN (fichnom, NF_NOWRITE,nid)
256	IF (ierr.NE.NF_NOERR) THEN
257	write(6,*)' Problem opening file:',fichnom
258	write(6,*)' ierr = ', ierr
259	CALL ABORT
260	ENDIF
261	tab0 = 50
262	Lmodif = 1
263
264	c OR open start and startfi files
265	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
266	else
267	write(,) 'Creating start files from:'
268	write(,) './start.nc and ./startfi.nc'
269	write(,)
270	fichnom = 'start.nc'
271	ierr = NF_OPEN (fichnom, NF_NOWRITE,nid_dyn)
272	IF (ierr.NE.NF_NOERR) THEN
273	write(6,*)' Problem opening file:',fichnom
274	write(6,*)' ierr = ', ierr
275	CALL ABORT
276	ENDIF
277
278	fichnom = 'startfi.nc'
279	ierr = NF_OPEN (fichnom, NF_NOWRITE,nid_fi)
280	IF (ierr.NE.NF_NOERR) THEN
281	write(6,*)' Problem opening file:',fichnom
282	write(6,*)' ierr = ', ierr
283	CALL ABORT
284	ENDIF
285
286	tab0 = 0
287	Lmodif = 0
288
289	endif
290
291	c-----------------------------------------------------------------------
292	c Lecture du tableau des parametres du run (pour la dynamique)
293	c-----------------------------------------------------------------------
294
295	if (choix_1.eq.0) then
296
297	write(,) 'reading tab_cntrl START_ARCHIVE'
298	c
299	ierr = NF_INQ_VARID (nid, "controle", nvarid)
300	#ifdef NC_DOUBLE
301	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tab_cntrl)
302	#else
303	ierr = NF_GET_VAR_REAL(nid, nvarid, tab_cntrl)
304	#endif
305	c
306	else if (choix_1.eq.1) then
307
308	write(,) 'reading tab_cntrl START'
309	c
310	ierr = NF_INQ_VARID (nid_dyn, "controle", nvarid)
311	#ifdef NC_DOUBLE
312	ierr = NF_GET_VAR_DOUBLE(nid_dyn, nvarid, tab_cntrl)
313	#else
314	ierr = NF_GET_VAR_REAL(nid_dyn, nvarid, tab_cntrl)
315	#endif
316	c
317	write(,) 'reading tab_cntrl STARTFI'
318	c
319	ierr = NF_INQ_VARID (nid_fi, "controle", nvarid)
320	#ifdef NC_DOUBLE
321	ierr = NF_GET_VAR_DOUBLE(nid_fi, nvarid, tab_cntrl_bis)
322	#else
323	ierr = NF_GET_VAR_REAL(nid_fi, nvarid, tab_cntrl_bis)
324	#endif
325	c
326	do i=1,50
327	tab_cntrl(i+50)=tab_cntrl_bis(i)
328	enddo
329	write(,) 'printing tab_cntrl', tab_cntrl
330	do i=1,100
331	write(,) i,tab_cntrl(i)
332	enddo
333
334	endif
335	c-----------------------------------------------------------------------
336	c Initialisation des constantes dynamique
337	c-----------------------------------------------------------------------
338
339	kappa = tab_cntrl(9)
340	etot0 = tab_cntrl(12)
341	ptot0 = tab_cntrl(13)
342	ztot0 = tab_cntrl(14)
343	stot0 = tab_cntrl(15)
344	ang0 = tab_cntrl(16)
345	write(,) "Newstart: kappa,etot0,ptot0,ztot0,stot0,ang0"
346	write(,) kappa,etot0,ptot0,ztot0,stot0,ang0
347
348	c-----------------------------------------------------------------------
349	c Lecture du tab_cntrl et initialisation des constantes physiques
350	c - pour start: Lmodif = 0 => pas de modifications possibles
351	c (modif dans le tabfi de readfi + loin)
352	c - pour start_archive: Lmodif = 1 => modifications possibles
353	c-----------------------------------------------------------------------
354	if (choix_1.eq.0) then
355	! tabfi requires that input file be first opened by open_startphy(fichnom)
356	fichnom = 'start_archive.nc'
357	call open_startphy(fichnom)
358	call tabfi (nid,Lmodif,tab0,day_ini,lllm,p_rad,
359	. p_omeg,p_g,p_mugaz,p_daysec,time)
360	else if (choix_1.eq.1) then
361	fichnom = 'startfi.nc'
362	call open_startphy(fichnom)
363	call tabfi (nid_fi,Lmodif,tab0,day_ini,lllm,p_rad,
364	. p_omeg,p_g,p_mugaz,p_daysec,time)
365	endif
366
367	rad = p_rad
368	omeg = p_omeg
369	g = p_g
370	mugaz = p_mugaz
371	daysec = p_daysec
372
373
374	c=======================================================================
375	c INITIALISATIONS DIVERSES
376	c=======================================================================
377
378	day_step=180 !?! Note: day_step is a common in "control.h"
379	CALL defrun_new( 99, .TRUE. )
380	dtvr = daysec/REAL(day_step)
381	CALL iniconst
382	CALL inigeom
383	idum=-1
384	idum=0
385
386	! Initialize the physics
387	CALL iniphysiq(iim,jjm,llm,
388	& (jjm-1)*iim+2,comm_lmdz,
389	& daysec,day_ini,dtphys,
390	& rlatu,rlatv,rlonu,rlonv,
391	& aire,cu,cv,rad,g,r,cpp,
392	& 1)
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	relief="mola"
407	c enddo
408
409	CALL datareadnc(relief,phis,alb,surfith,z0S,
410	& zmeaS,zstdS,zsigS,zgamS,ztheS,
411	& hmonsS,summitS,baseS,zavgS)
412
413	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
414	! CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ith,ithfi)
415	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,surfith,surfithfi)
416	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
417	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,z0S,z0)
418	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zmeaS,zmea)
419	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zstdS,zstd)
420	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zsigS,zsig)
421	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zgamS,zgam)
422	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ztheS,zthe)
423	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,hmonsS,hmons)
424	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,summitS,summit)
425	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,baseS,base)
426	! CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zavgS,zavg)
427
428	endif ! of if (choix_1.ne.1)
429
430
431	c=======================================================================
432	c Lecture des fichiers (start ou start_archive)
433	c=======================================================================
434
435	if (choix_1.eq.0) then
436
437	write(,) 'Reading file START_ARCHIVE'
438	CALL lect_start_archive(ngridmx,llm,nqtot,
439	& date,tsurf,tsoil,albedo,emis,q2,
440	& t,ucov,vcov,ps,teta,phisold_newgrid,q,qsurf,
441	& tauscaling,totcloudfrac,surfith,nid,watercap)
442	write(,) "OK, read start_archive file"
443	! copy soil thermal inertia
444	ithfi(:,:)=inertiedat(:,:)
445
446	ierr= NF_CLOSE(nid)
447
448	else if (choix_1.eq.1) then ! c'est l'appel a tabfi de phyeta0 qui
449	! permet de changer les valeurs du
450	! tab_cntrl Lmodif=1
451	tab0=0
452	Lmodif=1 ! Lmodif set to 1 to allow modifications in phyeta0
453	write(,) 'Reading file START'
454	fichnom = 'start.nc'
455	CALL dynetat0(fichnom,vcov,ucov,teta,q,masse,
456	& ps,phis,time)
457
458	write(,) 'Reading file STARTFI'
459	fichnom = 'startfi.nc'
460	CALL phyetat0 (fichnom,tab0,Lmodif,nsoilmx,ngridmx,llm,nqtot,
461	& day_ini,time,tsurf,tsoil,albedo,emis,
462	& q2,qsurf,tauscaling,totcloudfrac,
463	& wstar,watercap,def_slope,def_slope_mean,subslope_dist)
464
465	! copy albedo and soil thermal inertia
466	do i=1,ngridmx
467	albfi(i) = albedodat(i)
468	do j=1,nsoilmx
469	ithfi(i,j) = inertiedat(i,j)
470	enddo
471	! build a surfithfi(:) using 1st layer of ithfi(:), which might
472	! be neede later on if reinitializing soil thermal inertia
473	surfithfi(i)=ithfi(i,1)
474	enddo
475
476	else
477	CALL exit(1)
478	endif
479
480	dtvr = daysec/REAL(day_step)
481	dtphys = dtvr * REAL(iphysiq)
482
483	c=======================================================================
484	c
485	c=======================================================================
486	! If tracer names follow 'old' convention (q01, q02, ...) then
487	! rename them
488	count=0
489	do iq=1,nqtot
490	txt=" "
491	write(txt,'(a1,i2.2)') 'q',iq
492	if (txt.eq.tname(iq)) then
493	count=count+1
494	endif
495	enddo ! of do iq=1,nqtot
496
497	! initialize tracer names noms(:) and indexes (igcm_co2, igcm_h2o_vap, ...)
498	call initracer(ngridmx,nqtot,qsurf)
499
500	if (count.eq.nqtot) then
501	write(,) 'Newstart: updating tracer names'
502	! copy noms(:) to tname(:) to have matching tracer names in physics
503	! and dynamics
504	tname(1:nqtot)=noms(1:nqtot)
505	endif
506
507	c=======================================================================
508	c
509	c=======================================================================
510
511	do ! infinite loop on list of changes
512
513	write(,)
514	write(,)
515	write(,) 'List of possible changes :'
516	write(,) '~~~~~~~~~~~~~~~~~~~~~~~~~~'
517	write(,)
518	write(,) 'flat : no topography ("aquaplanet")'
519	write(,) 'bilball : uniform albedo and thermal inertia'
520	write(,) 'z0 : set a uniform surface roughness length'
521	write(,) 'coldspole : cold subsurface and high albedo at
522	$ S.Pole'
523	write(,) 'qname : change tracer name'
524	write(,) 'q=0 : ALL tracer =zero'
525	write(,) 'q=factor : change tracer value by a multiplicative
526	& factor'
527	write(,) 'q=x : give a specific uniform value to one
528	$ tracer'
529	write(,) 'q=profile : specify a profile for a tracer'
530	write(,) 'freedust : rescale dust to a true value'
531	write(,) 'ini_q : tracers initialization for chemistry
532	$ and water vapour'
533	write(,) 'ini_q-h2o : tracers initialization for chemistry
534	$ only'
535	write(,) 'composition : change atm main composition: CO2,N2,Ar,
536	$ O2,CO'
537	write(,) 'inihdo : initialize HDO'
538	write(,) 'ini_h2osurf : reinitialize surface water ice '
539	write(,) 'noglacier : Remove tropical H2O ice if \|lat\|<45'
540	write(,) 'watercapn : H20 ice on permanent N polar cap '
541	write(,) 'watercaps : H20 ice on permanent S polar cap '
542	write(,) 'wetstart : start with a wet atmosphere'
543	write(,) 'isotherm : Isothermal Temperatures, wind set to
544	$ zero'
545	write(,) 'co2ice=0 : remove CO2 polar cap i.e.
546	$ qsurf(co2)=0 '
547	write(,) 'ptot : change total pressure'
548	write(,) 'therm_ini_s : set soil thermal inertia to reference
549	$ surface values'
550	write(,) 'subsoilice_n : put deep underground ice layer in
551	$ northern hemisphere'
552	write(,) 'subsoilice_s : put deep underground ice layer in
553	$ southern hemisphere'
554	write(,) 'mons_ice : put underground ice layer according
555	$ to MONS derived data'
556	write(,) 'nslope : Change the number of subgrid scale
557	$ slope'
558
559	write(,)
560	write(,) 'Change to perform ?'
561	write(,) ' (enter keyword or return to end)'
562	write(,)
563
564	read(*,fmt='(a20)') modif
565	if (modif(1:1) .eq. ' ') exit ! exit loop on changes
566
567	write(,)
568	write(,) trim(modif) , ' : '
569
570	c 'flat : no topography ("aquaplanet")'
571	c -------------------------------------
572	if (trim(modif) .eq. 'flat') then
573	c set topo to zero
574	phis(:,:)=0
575	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
576	write(,) 'topography set to zero.'
577	write(,) 'WARNING : the subgrid topography parameters',
578	& ' were not set to zero ! => set calllott to F'
579
580	c Choice for surface pressure
581	yes=' '
582	do while ((yes.ne.'y').and.(yes.ne.'n'))
583	write(,) 'Do you wish to choose homogeneous surface',
584	& 'pressure (y) or let newstart interpolate ',
585	& ' the previous field (n)?'
586	read(*,fmt='(a)') yes
587	end do
588	if (yes.eq.'y') then
589	flagps0=.true.
590	write(,) 'New value for ps (Pa) ?'
591	201 read(,,iostat=ierr) patm
592	if(ierr.ne.0) goto 201
593	write(,)
594	write(,) ' new ps everywhere (Pa) = ', patm
595	write(,)
596	do j=1,jjp1
597	do i=1,iip1
598	ps(i,j)=patm
599	enddo
600	enddo
601	end if
602
603	c bilball : albedo, inertie thermique uniforme
604	c --------------------------------------------
605	else if (trim(modif) .eq. 'bilball') then
606	write(,) 'constante albedo and iner.therm:'
607	write(,) 'New value for albedo (ex: 0.25) ?'
608	101 read(,,iostat=ierr) alb_bb
609	if(ierr.ne.0) goto 101
610	write(,)
611	write(,) ' uniform albedo (new value):',alb_bb
612	write(,)
613
614	write(,) 'New value for thermal inertia (eg: 247) ?'
615	102 read(,,iostat=ierr) ith_bb
616	if(ierr.ne.0) goto 102
617	write(,) 'uniform thermal inertia (new value):',ith_bb
618	DO j=1,jjp1
619	DO i=1,iip1
620	alb(i,j) = alb_bb ! albedo
621	do isoil=1,nsoilmx
622	ith(i,j,isoil) = ith_bb ! thermal inertia
623	enddo
624	END DO
625	END DO
626	! CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ith,ithfi)
627	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
628	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
629
630	! also reset surface roughness length to default value
631	write(,) 'surface roughness length set to:',z0_default,' m'
632	z0(:)=z0_default
633
634	! z0 : set surface roughness length to a constant value
635	! -----------------------------------------------------
636	else if (trim(modif) .eq. 'z0') then
637	write(,) 'set a uniform surface roughness length'
638	write(,) ' value for z0_default (ex: ',z0_default,')?'
639	ierr=1
640	do while (ierr.ne.0)
641	read(,,iostat=ierr) z0_default
642	enddo
643	z0(:)=z0_default
644
645	c coldspole : sous-sol de la calotte sud toujours froid
646	c -----------------------------------------------------
647	else if (trim(modif) .eq. 'coldspole') then
648	write(,)'new value for the subsurface temperature',
649	& ' beneath the permanent southern polar cap ? (eg: 141 K)'
650	103 read(,,iostat=ierr) tsud
651	if(ierr.ne.0) goto 103
652	write(,)
653	write(,) ' new value of the subsurface temperature:',tsud
654	c nouvelle temperature sous la calotte permanente
655	do islope=1,nslope
656	do l=2,nsoilmx
657	tsoil(ngridmx,l,islope) = tsud
658	end do
659	enddo
660
661
662	write(,)'new value for the albedo',
663	& 'of the permanent southern polar cap ? (eg: 0.75)'
664	104 read(,,iostat=ierr) albsud
665	if(ierr.ne.0) goto 104
666	write(,)
667
668	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
669	c Option 1: only the albedo of the pole is modified :
670	albfi(ngridmx)=albsud
671	write(,) 'ig=',ngridmx,' albedo perennial cap ',
672	& albfi(ngridmx)
673
674	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
675	c Option 2 A haute resolution : coordonnee de la vrai calotte ~
676	c DO j=1,jjp1
677	c DO i=1,iip1
678	c ig=1+(j-2)*iim +i
679	c if(j.eq.1) ig=1
680	c if(j.eq.jjp1) ig=ngridmx
681	c if ((rlatu(j)*180./pi.lt.-84.).and.
682	c & (rlatu(j)*180./pi.gt.-91.).and.
683	c & (rlonv(i)*180./pi.gt.-91.).and.
684	c & (rlonv(i)*180./pi.lt.0.)) then
685	cc albedo de la calotte permanente fixe a albsud
686	c alb(i,j)=albsud
687	c write(,) 'lat=',rlatu(j)*180./pi,
688	c & ' lon=',rlonv(i)*180./pi
689	cc fin de la condition sur les limites de la calotte permanente
690	c end if
691	c ENDDO
692	c ENDDO
693	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
694
695	c CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
696
697
698	c ptot : Modification of the total pressure: ice + current atmosphere
699	c -------------------------------------------------------------------
700	else if (trim(modif) .eq. 'ptot') then
701
702	c calcul de la pression totale glace + atm actuelle
703	patm=0.
704	airetot=0.
705	pcap=0.
706	DO j=1,jjp1
707	DO i=1,iim
708	ig=1+(j-2)*iim +i
709	if(j.eq.1) ig=1
710	if(j.eq.jjp1) ig=ngridmx
711	patm = patm + ps(i,j)*aire(i,j)
712	airetot= airetot + aire(i,j)
713	DO islope=1,nslope
714	pcap = pcap + aire(i,j)qsurf(ig,igcm_co2,islope)g
715	& subslope_dist(ig,islope)/cos(def_slope_mean(islope)pi/180.)
716	ENDDO
717	ENDDO
718	ENDDO
719	ptoto = pcap + patm
720
721	print*,'Current total pressure at surface (co2 ice + atm) ',
722	& ptoto/airetot
723
724	print*,'new value?'
725	read(,) ptotn
726	ptotn=ptotn*airetot
727	patmn=ptotn-pcap
728	print*,'ptoto,patm,ptotn,patmn'
729	print*,ptoto,patm,ptotn,patmn
730	print*,'Mult. factor for pressure (atm only)', patmn/patm
731	do j=1,jjp1
732	do i=1,iip1
733	ps(i,j)=ps(i,j)*patmn/patm
734	enddo
735	enddo
736
737	c Correction pour la conservation des traceurs
738	yes=' '
739	do while ((yes.ne.'y').and.(yes.ne.'n'))
740	write(,) 'Do you wish to conserve tracer total mass (y)',
741	& ' or tracer mixing ratio (n) ?'
742	read(*,fmt='(a)') yes
743	end do
744
745	if (yes.eq.'y') then
746	write(,) 'OK : conservation of tracer total mass'
747	DO iq =1, nqtot
748	DO l=1,llm
749	DO j=1,jjp1
750	DO i=1,iip1
751	q(i,j,l,iq)=q(i,j,l,iq)*patm/patmn
752	ENDDO
753	ENDDO
754	ENDDO
755	ENDDO
756	else
757	write(,) 'OK : conservation of tracer mixing ratio'
758	end if
759
760	c qname : change tracer name
761	c --------------------------
762	else if (trim(modif).eq.'qname') then
763	yes='y'
764	do while (yes.eq.'y')
765	write(,) 'Which tracer name do you want to change ?'
766	do iq=1,nqtot
767	write(*,'(i3,a3,a20)')iq,' : ',trim(tname(iq))
768	enddo
769	write(*,'(a35,i3)')
770	& '(enter tracer number; between 1 and ',nqtot
771	write(,)' or any other value to quit this option)'
772	read(,) iq
773	if ((iq.ge.1).and.(iq.le.nqtot)) then
774	write(,)'Change tracer name ',trim(tname(iq)),' to ?'
775	read(,) txt
776	tname(iq)=txt
777	write(,)'Do you want to change another tracer name (y/n)?'
778	read(*,'(a)') yes
779	else
780	! inapropiate value of iq; quit this option
781	yes='n'
782	endif ! of if ((iq.ge.1).and.(iq.le.nqtot))
783	enddo ! of do while (yes.ne.'y')
784
785	c q=0 : set tracers to zero
786	c -------------------------
787	else if (trim(modif) .eq. 'q=0') then
788	c mise a 0 des q (traceurs)
789	write(,) 'Tracers set to 0 (1.E-30 in fact)'
790	DO iq =1, nqtot
791	DO l=1,llm
792	DO j=1,jjp1
793	DO i=1,iip1
794	q(i,j,l,iq)=1.e-30
795	ENDDO
796	ENDDO
797	ENDDO
798	ENDDO
799
800	c set surface tracers to zero
801	DO iq =1, nqtot
802	DO ig=1,ngridmx
803	DO islope=1,nslope
804	qsurf(ig,iq,islope)=0.
805	ENDDO
806	ENDDO
807	ENDDO
808
809	c q=factor : change value of tracer by a multiplicative factor
810	c ------------------------------------------------------------
811	else if (trim(modif) .eq. 'q=factor') then
812	write(,) 'Which tracer do you want to modify ?'
813	do iq=1,nqtot
814	write(,)iq,' : ',trim(tname(iq))
815	enddo
816	write(,) '(choose between 1 and ',nqtot,')'
817	read(,) iq
818	if ((iq.lt.1).or.(iq.gt.nqtot)) then
819	! wrong value for iq, go back to menu
820	write(,) "wrong input value:",iq
821	cycle
822	endif
823	write(,)"factor to multiply current mixing ratio by?"
824	read(,) val
825
826	q(1:iip1,1:jjp1,1:llm,iq)=q(1:iip1,1:jjp1,1:llm,iq)*val
827	qsurf(1:ngridmx,iq,:)=qsurf(1:ngridmx,iq,:)*val
828
829	c q=x : initialise tracer manually
830	c --------------------------------
831	else if (trim(modif) .eq. 'q=x') then
832	write(,) 'Which tracer do you want to modify ?'
833	do iq=1,nqtot
834	write(,)iq,' : ',trim(tname(iq))
835	enddo
836	write(,) '(choose between 1 and ',nqtot,')'
837	read(,) iq
838	if ((iq.lt.1).or.(iq.gt.nqtot)) then
839	! wrong value for iq, go back to menu
840	write(,) "wrong input value:",iq
841	cycle
842	endif
843	write(,)'mixing ratio of tracer ',trim(tname(iq)),
844	& ' ? (kg/kg)'
845	read(,) val
846	DO l=1,llm
847	DO j=1,jjp1
848	DO i=1,iip1
849	q(i,j,l,iq)=val
850	ENDDO
851	ENDDO
852	ENDDO
853	write(,) 'SURFACE value of tracer ',trim(tname(iq)),
854	& ' ? (kg/m2)'
855	read(,) val
856	DO ig=1,ngridmx
857	DO islope=1,nslope
858	qsurf(ig,iq,islope)=val
859	ENDDO
860	ENDDO
861
862	c q=profile : initialize tracer with a given profile
863	c --------------------------------------------------
864	else if (trim(modif) .eq. 'q=profile') then
865	write(,) 'Tracer profile will be sought in ASCII file'
866	write(,) "'profile_tracer' where 'tracer' is tracer name"
867	write(,) "(one value per line in file; starting with"
868	write(,) "surface value, the 1st atmospheric layer"
869	write(,) "followed by 2nd, etc. up to top of atmosphere)"
870	write(,) 'Which tracer do you want to set?'
871	do iq=1,nqtot
872	write(,)iq,' : ',trim(tname(iq))
873	enddo
874	write(,) '(choose between 1 and ',nqtot,')'
875	read(,) iq
876	if ((iq.lt.1).or.(iq.gt.nqtot)) then
877	! wrong value for iq, go back to menu
878	write(,) "wrong input value:",iq
879	cycle
880	endif
881	! look for input file 'profile_tracer'
882	txt="profile_"//trim(tname(iq))
883	open(41,file=trim(txt),status='old',form='formatted',
884	& iostat=ierr)
885	if (ierr.eq.0) then
886	! OK, found file 'profile_...', load the profile
887	do l=1,llm+1
888	read(41,*,iostat=ierr) profile(l)
889	if (ierr.ne.0) then ! something went wrong
890	exit ! quit loop
891	endif
892	enddo
893	if (ierr.eq.0) then
894	! initialize tracer values
895	do islope=1,nslope
896	qsurf(:,iq,islope)=profile(1)
897	enddo
898	do l=1,llm
899	q(:,:,l,iq)=profile(l+1)
900	enddo
901	write(,)'OK, tracer ',trim(tname(iq)),
902	& ' initialized ','using values from file ',trim(txt)
903	else
904	write(,)'problem reading file ',trim(txt),' !'
905	write(,)'No modifications to tracer ',trim(tname(iq))
906	endif
907	else
908	write(,)'Could not find file ',trim(txt),' !'
909	write(,)'No modifications to tracer ',trim(tname(iq))
910	endif
911
912	c convert dust from virtual to true values
913	c --------------------------------------------------
914	else if (trim(modif) .eq. 'freedust') then
915	if (minval(tauscaling) .lt. 0) then
916	write(,) 'WARNING conversion factor negative'
917	write(,) 'This is probably because it was not present
918	&in the file'
919	write(,) 'A constant conversion is used instead.'
920	tauscaling(:) = 1.e-3
921	endif
922	CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,tauscaling,tauscadyn)
923	do l=1,llm
924	do j=1,jjp1
925	do i=1,iip1
926	if (igcm_dust_number .ne. 0)
927	& q(i,j,l,igcm_dust_number) =
928	& q(i,j,l,igcm_dust_number) * tauscadyn(i,j)
929	if (igcm_dust_mass .ne. 0)
930	& q(i,j,l,igcm_dust_mass) =
931	& q(i,j,l,igcm_dust_mass) * tauscadyn(i,j)
932	if (igcm_ccn_number .ne. 0)
933	& q(i,j,l,igcm_ccn_number) =
934	& q(i,j,l,igcm_ccn_number) * tauscadyn(i,j)
935	if (igcm_ccn_mass .ne. 0)
936	& q(i,j,l,igcm_ccn_mass) =
937	& q(i,j,l,igcm_ccn_mass) * tauscadyn(i,j)
938	end do
939	end do
940	end do
941
942	tauscaling(:) = 1.
943
944	! We want to have the very same value at lon -180 and lon 180
945	do l = 1,llm
946	do j = 1,jjp1
947	do iq = 1,nqtot
948	q(iip1,j,l,iq) = q(1,j,l,iq)
949	end do
950	end do
951	end do
952
953	write(,) 'done rescaling to true vale'
954
955	c ini_q : Initialize tracers for chemistry
956	c -----------------------------------------------
957	else if (trim(modif) .eq. 'ini_q') then
958	flagh2o = 1
959	flagthermo = 0
960	yes=' '
961	c For more than 32 layers, possible to initiate thermosphere only
962	if (llm.gt.32) then
963	do while ((yes.ne.'y').and.(yes.ne.'n'))
964	write(,)'',
965	& 'initialisation for thermosphere only? (y/n)'
966	read(*,fmt='(a)') yes
967	if (yes.eq.'y') then
968	flagthermo=1
969	else
970	flagthermo=0
971	endif
972	enddo
973	endif
974
975	call inichim_newstart(ngridmx, nqtot, q, qsurf, ps,
976	& flagh2o, flagthermo)
977
978	! We want to have the very same value at lon -180 and lon 180
979	do l = 1,llm
980	do j = 1,jjp1
981	do iq = 1,nqtot
982	q(iip1,j,l,iq) = q(1,j,l,iq)
983	end do
984	end do
985	end do
986
987	write(,) 'inichim_newstart: chemical species and
988	$ water vapour initialised'
989
990	c ini_q-h2o : as above except for the water vapour tracer
991	c ------------------------------------------------------
992	else if (trim(modif) .eq. 'ini_q-h2o') then
993	flagh2o = 0
994	flagthermo = 0
995	yes=' '
996	! for more than 32 layers, possible to initiate thermosphere only
997	if(llm.gt.32) then
998	do while ((yes.ne.'y').and.(yes.ne.'n'))
999	write(,)'',
1000	& 'initialisation for thermosphere only? (y/n)'
1001	read(*,fmt='(a)') yes
1002	if (yes.eq.'y') then
1003	flagthermo=1
1004	else
1005	flagthermo=0
1006	endif
1007	enddo
1008	endif
1009
1010	call inichim_newstart(ngridmx, nqtot, q, qsurf, ps,
1011	& flagh2o, flagthermo)
1012
1013	! We want to have the very same value at lon -180 and lon 180
1014	do l = 1,llm
1015	do j = 1,jjp1
1016	do iq = 1,nqtot
1017	q(iip1,j,l,iq) = q(1,j,l,iq)
1018	end do
1019	end do
1020	end do
1021
1022	write(,) 'inichim_newstart: chemical species initialised
1023	$ (except water vapour)'
1024
1025	c inihdo : initialize HDO with user D/H value
1026	c --------------------------------------------------------
1027	else if (trim(modif) .eq. 'inihdo') then
1028	! check that there is indeed a water vapor tracer
1029	if (igcm_h2o_vap.eq.0) then
1030	write(,) "No water vapour tracer! Can't use this option"
1031	stop
1032	endif
1033
1034	write(,)'Input D/H ratio (in SMOW)'
1035	write(,)'If value is <0 then HDO=H2O'
1036	303 read(,, iostat=ierr) DoverH
1037	if(ierr.ne.0) goto 303
1038
1039	DoverH = DoverH * 2 * 155.76e-6
1040	if (DoverH.lt.0.0) then
1041	DoverH = 1.
1042	endif
1043	!D/H (SMOW) = 155.76e-6 so HDO/H2O is twice that
1044
1045	do ig=1,ngridmx
1046	do islope=1,nslope
1047	qsurf(ig,igcm_h2o_ice,islope)=
1048	& max(0.,qsurf(ig,igcm_h2o_ice,islope))
1049	enddo
1050	end do
1051
1052	! Update the hdo tracers
1053	q(1:iip1,1:jjp1,1:llm,igcm_hdo_vap)
1054	& =q(1:iip1,1:jjp1,1:llm,igcm_h2o_vap)* DoverH
1055	q(1:iip1,1:jjp1,1:llm,igcm_hdo_ice)
1056	& =q(1:iip1,1:jjp1,1:llm,igcm_h2o_ice)* DoverH
1057
1058	qsurf(1:ngridmx,igcm_hdo_ice,:)
1059	& =qsurf(1:ngridmx,igcm_h2o_ice,:)*DoverH
1060
1061
1062
1063	c composition : change main composition: CO2,N2,Ar,O2,CO (FF 03/2014)
1064	c --------------------------------------------------------
1065	else if (trim(modif) .eq. 'composition') then
1066	write(,) "Lat (degN) lon (degE) of the reference site ?"
1067	write(,) "e.g. MSL : -4.5 137. "
1068	301 read(,,iostat=ierr) latref, lonref
1069	if(ierr.ne.0) goto 301
1070
1071
1072	! Select GCM point close to reference site
1073	dlonmin =90.
1074	DO i=1,iip1-1
1075	if (abs(rlonv(i)*180./pi -lonref).lt.dlonmin)then
1076	iref=i
1077	dlonmin=abs(rlonv(i)*180./pi -lonref)
1078	end if
1079	ENDDO
1080	dlatmin =45.
1081	DO j=1,jjp1
1082	if (abs(rlatu(j)*180./pi -latref).lt.dlatmin)then
1083	jref=j
1084	dlatmin=abs(rlatu(j)*180./pi -latref)
1085	end if
1086	ENDDO
1087	write(,) "In GCM : lat= " , rlatu(jref)*180./pi
1088	write(,) "In GCM : lon= " , rlonv(iref)*180./pi
1089	write(,)
1090
1091	! Compute air molar mass at reference site
1092	Smmr=0.
1093	Sn = 0.
1094	write(,) 'igcm_co2 = ', igcm_co2
1095	write(,) 'igcm_n2 = ', igcm_n2
1096	write(,) 'igcm_ar = ', igcm_ar
1097	write(,) 'igcm_o2 = ', igcm_o2
1098	write(,) 'igcm_co = ', igcm_co
1099	write(,) noms
1100	do iq=1,nqtot
1101	if ((iq.eq.igcm_co2).or.(iq.eq.igcm_n2)
1102	& .or. (iq.eq.igcm_ar).or.(iq.eq.igcm_o2)
1103	& .or. (iq.eq.igcm_co)) then
1104	Smmr=Smmr+q(iref,jref,1,iq)
1105	Sn=Sn+q(iref,jref,1,iq)/mmol(iq)
1106	end if
1107	end do
1108	! Special case : implicit non-co2 gases ! JN 11/2019
1109	if ((igcm_n2.eq.0) .or. (igcm_ar.eq.0)) then
1110	write(,) "Warning : non-co2 gases are implicit : "
1111	write(,) "At reference site : "
1112	! write(,) "q= ", q(iref, jref, 1,igcm_co2)
1113	write(,) "Sum of mass mix. ratio (ie MMR(co2))=",Smmr
1114	Mair_old = 44.0Smmr + 33.87226017157708(1-Smmr)
1115
1116	! 33.87226017157708 is the
1117	! molar mass of non-co2 atmosphere measured by MSL at Ls ~184
1118
1119	else
1120	! Assume co2/n2/ar/o2/co are available
1121	Mair_old=(q(iref,jref,1,igcm_co2)*mmol(igcm_co2)
1122	& +q(iref,jref,1,igcm_n2)*mmol(igcm_n2)
1123	& +q(iref,jref,1,igcm_ar)*mmol(igcm_ar)
1124	& +q(iref,jref,1,igcm_o2)*mmol(igcm_o2)
1125	& +q(iref,jref,1,igcm_co)*mmol(igcm_co))/Smmr
1126	end if
1127
1128	write(,)
1129	& "Air molar mass (g/mol) at reference site= ",Mair_old
1130
1131	! Ask for new volume mixing ratio at reference site
1132	Svmr =0.
1133	Sn =0.
1134	coefvmr(igcm_co2)=1.
1135
1136	do iq=1,nqtot
1137	coefvmr(iq) = 1.
1138	if ((iq.eq.igcm_n2).or.(iq.eq.igcm_ar)
1139	& .or. (iq.eq.igcm_o2).or.(iq.eq.igcm_co)) then
1140
1141	vmr_old=q(iref,jref,1,iq)*Mair_old/mmol(iq)
1142	write(,) "Previous vmr("//trim(tname(iq))//")= ", vmr_old
1143
1144	if (iq.eq.igcm_n2) then
1145	write(,) "New vmr(n2)? (MSL: 2.8e-02 at Ls~180,",
1146	& " Trainer et al. 2019)"
1147	endif
1148	if (iq.eq.igcm_ar) then
1149	write(,) "New vmr(ar)? (MSL: 2.1e-02 at Ls~180)"
1150	endif
1151	if (iq.eq.igcm_o2) then
1152	write(,) "New vmr(o2)? (MSL: 1.7e-03 at Ls~180)"
1153	endif
1154	if (iq.eq.igcm_co) then
1155	write(,) "New vmr(co)? (ACS: 1e-03 at Ls~180)"
1156	endif
1157	302 read(,,iostat=ierr) vmr_new
1158	if(ierr.ne.0) goto 302
1159	write(,) "New vmr("//trim(tname(iq))//")= ",vmr_new
1160	write(,)
1161	coefvmr(iq) = vmr_new/vmr_old
1162	Svmr=Svmr+vmr_new
1163	Sn=Sn+vmr_new*mmol(iq)
1164	end if
1165	enddo ! of do iq=1,nqtot
1166
1167	! Special case : implicit non-co2 gases JN 11/2019
1168	if ((igcm_n2.eq.0) .or. (igcm_ar.eq.0)) then
1169	write(,) "Warning : non-co2 gases are implicit"
1170	vmr_old=q(iref,jref,1,igcm_co2)*Mair_old/mmol(igcm_co2)
1171	write(,) "Previous vmr(co2)=", vmr_old
1172	write(,) "New vmr(co2) ? (MSL: 0.947 at Ls~180)",
1173	& " Trainer et al. 2019)"
1174	666 read(,,iostat=ierr) vmr_new
1175	if(ierr.ne.0) goto 666
1176	coefvmr(igcm_co2) = vmr_new/vmr_old
1177	Svmr=Svmr+vmr_new
1178	Sn=vmr_new*mmol(igcm_co2) + (1-vmr_new)
1179	& *33.87226017157708 ! Molar mass of non-co2 atm from MSL
1180	end if
1181	! Estimation of new Air molar mass at reference site (assuming vmr_co2 = 1-Svmr)
1182	Mair_new = Sn + (1-Svmr)*mmol(igcm_co2)
1183	! Estimation of new Air molar mass when non-co2 gases are implicit
1184	if ((igcm_n2.eq.0) .or. (igcm_ar.eq.0)) then
1185	Mair_new=vmr_new*mmol(igcm_co2) + (1-vmr_new)
1186	& *33.87226017157708 ! Molar mass of non-co2 atm from MSL
1187	write(,)
1188	& "We consider non-co2 gases vmr measured from Curiosity"
1189	end if
1190	write(,)
1191	& "NEW Air molar mass (g/mol) at reference site= ",Mair_new
1192
1193	! Compute mass mixing ratio changes
1194	do iq=1,nqtot
1195	if ((iq.eq.igcm_co2).or.(iq.eq.igcm_n2).or.(iq.eq.igcm_ar)
1196	& .or. (iq.eq.igcm_o2).or.(iq.eq.igcm_co)) then
1197	write(,) "Everywhere mmr("//trim(tname(iq))//
1198	& ") is multiplied by ",coefvmr(iq)*Mair_old/Mair_new
1199	end if
1200	end do
1201
1202	! Recompute mass mixing ratios everywhere, and adjust mmr of most abundant species
1203	! to keep sum of mmr constant.
1204	do l=1,llm
1205	do j=1,jjp1
1206	do i=1,iip1
1207	Smmr_old = 0.
1208	Smmr_new = 0.
1209	do iq=1,nqtot
1210	if ((iq.eq.igcm_co2).or.(iq.eq.igcm_n2)
1211	& .or.(iq.eq.igcm_ar)
1212	& .or. (iq.eq.igcm_o2).or.(iq.eq.igcm_co)
1213	& .or. (iq.eq.igcm_o) .or. (iq.eq. igcm_h2) ) then
1214	Smmr_old = Smmr_old + q(i,j,l,iq) ! sum of old mmr
1215	q(i,j,l,iq)=q(i,j,l,iq)coefvmr(iq)Mair_old/Mair_new
1216	Smmr_new = Smmr_new + q(i,j,l,iq) ! sum of new mmr
1217	end if
1218	enddo
1219	!iloc = maxloc(q(i,j,l,:))
1220	iqmax=0 ; maxq=0
1221	do iq=1,nqtot
1222	if ((iq.eq.igcm_co2).or.(iq.eq.igcm_n2)
1223	& .or.(iq.eq.igcm_ar)
1224	& .or. (iq.eq.igcm_o2).or.(iq.eq.igcm_co)
1225	& .or. (iq.eq.igcm_o) .or. (iq.eq. igcm_h2) ) then
1226	if (q(i,j,l,iq).gt.maxq) then
1227	maxq=q(i,j,l,iq)
1228	iqmax=iq
1229	endif
1230	endif
1231	enddo
1232	!iqmax = iloc(1)
1233	q(i,j,l,iqmax) = q(i,j,l,iqmax) + Smmr_old - Smmr_new
1234	enddo
1235	enddo
1236	enddo
1237
1238	write(,)
1239	& "The most abundant species is modified everywhere to keep "//
1240	& "sum of mmr constant"
1241	write(,) 'At reference site vmr(CO2)=',
1242	& q(iref,jref,1,igcm_co2)*Mair_new/mmol(igcm_co2)
1243	write(,) "Compared to MSL observation: vmr(CO2)= 0.947 "//
1244	& "at Ls=180"
1245
1246	Sn = q(iref,jref,1,igcm_co2)*Mair_new/mmol(igcm_co2)
1247	& + q(iref,jref,1,igcm_n2)*Mair_new/mmol(igcm_n2)
1248	& + q(iref,jref,1,igcm_ar)*Mair_new/mmol(igcm_ar)
1249	& + q(iref,jref,1,igcm_o2)*Mair_new/mmol(igcm_o2)
1250	& + q(iref,jref,1,igcm_co)*Mair_new/mmol(igcm_co)
1251
1252	write(,) 'Sum of volume mixing ratios = ', Sn
1253
1254	c wetstart : wet atmosphere with a north to south gradient
1255	c --------------------------------------------------------
1256	else if (trim(modif) .eq. 'wetstart') then
1257	! check that there is indeed a water vapor tracer
1258	if (igcm_h2o_vap.eq.0) then
1259	write(,) "No water vapour tracer! Can't use this option"
1260	stop
1261	endif
1262	DO l=1,llm
1263	DO j=1,jjp1
1264	DO i=1,iip1-1
1265	q(i,j,l,igcm_h2o_vap)=150.e-6 * (rlatu(j)+pi/2.) / pi
1266	ENDDO
1267	! We want to have the very same value at lon -180 and lon 180
1268	q(iip1,j,l,igcm_h2o_vap) = q(1,j,l,igcm_h2o_vap)
1269	ENDDO
1270	ENDDO
1271
1272	write(,) 'Water mass mixing ratio at north pole='
1273	* ,q(1,1,1,igcm_h2o_vap)
1274	write(,) '---------------------------south pole='
1275	* ,q(1,jjp1,1,igcm_h2o_vap)
1276
1277	c ini_h2osurf : reinitialize surface water ice
1278	c --------------------------------------------------
1279	else if (trim(modif) .eq. 'ini_h2osurf') then
1280	write(,)'max surface ice left?(e.g. 0.2 kg/m2=200microns)'
1281	207 read(,,iostat=ierr) val
1282	if(ierr.ne.0) goto 207
1283	write(,)'also set negative values of surf ice to 0'
1284	do ig=1,ngridmx
1285	do islope=1,nslope
1286	qsurf(ig,igcm_h2o_ice,islope)=
1287	& min(val,qsurf(ig,igcm_h2o_ice,islope))
1288	qsurf(ig,igcm_h2o_ice,islope)=
1289	& max(0.,qsurf(ig,igcm_h2o_ice,islope))
1290	enddo
1291	end do
1292
1293	c noglacier : remove tropical water ice (to initialize high res sim)
1294	c --------------------------------------------------
1295	else if (trim(modif) .eq. 'noglacier') then
1296	do ig=1,ngridmx
1297	j=(ig-2)/iim +2
1298	if(ig.eq.1) j=1
1299	write(,) 'OK: remove surface ice for \|lat\|<45'
1300	if (abs(rlatu(j)*180./pi).lt.45.) then
1301	do islope=1,nslope
1302	qsurf(ig,igcm_h2o_ice,islope)=0.
1303	enddo
1304	end if
1305	end do
1306
1307
1308	c watercapn : H20 ice on permanent northern cap
1309	c --------------------------------------------------
1310	else if (trim(modif) .eq. 'watercapn') then
1311	do ig=1,ngridmx
1312	j=(ig-2)/iim +2
1313	if(ig.eq.1) j=1
1314	if (rlatu(j)*180./pi.gt.80.) then
1315	do islope=1,nslope
1316	qsurf(ig,igcm_h2o_ice,islope)=1.e5
1317	write(,) 'ig=',ig,', islope=', islope,
1318	& ' H2O ice mass (kg/m2)= ',
1319	& qsurf(ig,igcm_h2o_ice,islope)
1320	write(,)' ==> Ice mesh South boundary (deg)= ',
1321	& rlatv(j)*180./pi
1322	enddo
1323	end if
1324	enddo
1325
1326	c watercaps : H20 ice on permanent southern cap
1327	c -------------------------------------------------
1328	else if (trim(modif) .eq. 'watercaps') then
1329	do ig=1,ngridmx
1330	j=(ig-2)/iim +2
1331	if(ig.eq.1) j=1
1332	if (rlatu(j)*180./pi.lt.-80.) then
1333	do islope=1,nslope
1334	qsurf(ig,igcm_h2o_ice,islope)=1.e5
1335	write(,) 'ig=',ig,', islope=', islope,
1336	& ' H2O ice mass (kg/m2)= ',
1337	& qsurf(ig,igcm_h2o_ice,islope)
1338	write(,)' ==> Ice mesh North boundary (deg)= ',
1339	& rlatv(j-1)*180./pi
1340	enddo
1341	end if
1342	enddo
1343
1344	c isotherm : Isothermal temperatures and no winds
1345	c ------------------------------------------------
1346	else if (trim(modif) .eq. 'isotherm') then
1347
1348	write(,)'Isothermal temperature of the atmosphere,
1349	& surface and subsurface'
1350	write(,) 'Value of this temperature ? :'
1351	203 read(,,iostat=ierr) Tiso
1352	if(ierr.ne.0) goto 203
1353
1354	do ig=1, ngridmx
1355	do islope=1,nslope
1356	tsurf(ig,islope) = Tiso
1357	enddo
1358	end do
1359	do l=2,nsoilmx
1360	do ig=1, ngridmx
1361	do islope=1,nslope
1362	tsoil(ig,l,islope) = Tiso
1363	enddo
1364	end do
1365	end do
1366	flagiso=.true.
1367	call initial0(llm*ip1jmp1,ucov)
1368	call initial0(llm*ip1jm,vcov)
1369	call initial0(ngridmx*(llm+1),q2)
1370
1371	c co2ice=0 : remove CO2 polar ice caps'
1372	c ------------------------------------------------
1373	else if (trim(modif) .eq. 'co2ice=0') then
1374	do ig=1,ngridmx
1375	do islope=1,nslope
1376	qsurf(ig,igcm_co2,islope)=0
1377	emis(ig,islope)=emis(ngridmx/2,islope)
1378	enddo
1379	end do
1380
1381	! therm_ini_s: (re)-set soil thermal inertia to reference surface values
1382	! ----------------------------------------------------------------------
1383
1384	else if (trim(modif).eq.'therm_ini_s') then
1385	! write(,)"surfithfi(1):",surfithfi(1)
1386	do isoil=1,nsoilmx
1387	inertiedat(1:ngridmx,isoil)=surfithfi(1:ngridmx)
1388	enddo
1389	do islope = 1,nslope
1390	inertiesoil(:,:,islope) = inertiedat(:,:)
1391	enddo
1392	write(,)'OK: Soil thermal inertia has been reset to referenc
1393	&e surface values'
1394	! write(,)"inertiedat(1,1):",inertiedat(1,1)
1395	ithfi(:,:)=inertiedat(:,:)
1396	! recast ithfi() onto ith()
1397	call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
1398	! Check:
1399	! do i=1,iip1
1400	! do j=1,jjp1
1401	! do isoil=1,nsoilmx
1402	! write(77,*) i,j,isoil," ",ith(i,j,isoil)
1403	! enddo
1404	! enddo
1405	! enddo
1406
1407	! subsoilice_n: Put deep ice layer in northern hemisphere soil
1408	! ------------------------------------------------------------
1409
1410	else if (trim(modif).eq.'subsoilice_n') then
1411
1412	write(,)'From which latitude (in deg.), up to the north pole,
1413	&should we put subterranean ice?'
1414	ierr=1
1415	do while (ierr.ne.0)
1416	read(,,iostat=ierr) val
1417	if (ierr.eq.0) then ! got a value
1418	! do a sanity check
1419	if((val.lt.0.).or.(val.gt.90)) then
1420	write(,)'Latitude should be between 0 and 90 deg. !!!'
1421	ierr=1
1422	else ! find corresponding jref (nearest latitude)
1423	! note: rlatu(:) contains decreasing values of latitude
1424	! starting from PI/2 to -PI/2
1425	do j=1,jjp1
1426	if ((rlatu(j)*180./pi.ge.val).and.
1427	& (rlatu(j+1)*180./pi.le.val)) then
1428	! find which grid point is nearest to val:
1429	if (abs(rlatu(j)*180./pi-val).le.
1430	& abs((rlatu(j+1)*180./pi-val))) then
1431	jref=j
1432	else
1433	jref=j+1
1434	endif
1435
1436	write(,)'Will use nearest grid latitude which is:',
1437	& rlatu(jref)*180./pi
1438	endif
1439	enddo ! of do j=1,jjp1
1440	endif ! of if((val.lt.0.).or.(val.gt.90))
1441	endif !of if (ierr.eq.0)
1442	enddo ! of do while
1443
1444	! Build layers() (as in soil_settings.F)
1445	val2=sqrt(mlayer(0)*mlayer(1))
1446	val3=mlayer(1)/mlayer(0)
1447	do isoil=1,nsoilmx
1448	layer(isoil)=val2(val3*(isoil-1))
1449	enddo
1450
1451	write(,)'At which depth (in m.) does the ice layer begin?'
1452	write(,)'(currently, the deepest soil layer extends down to:'
1453	& ,layer(nsoilmx),')'
1454	ierr=1
1455	do while (ierr.ne.0)
1456	read(,,iostat=ierr) val2
1457	! write(,)'val2:',val2,'ierr=',ierr
1458	if (ierr.eq.0) then ! got a value, but do a sanity check
1459	if(val2.gt.layer(nsoilmx)) then
1460	write(,)'Depth should be less than ',layer(nsoilmx)
1461	ierr=1
1462	endif
1463	if(val2.lt.layer(1)) then
1464	write(,)'Depth should be more than ',layer(1)
1465	ierr=1
1466	endif
1467	endif
1468	enddo ! of do while
1469
1470	! find the reference index iref the depth corresponds to
1471	! if (val2.lt.layer(1)) then
1472	! iref=1
1473	! else
1474	do isoil=1,nsoilmx-1
1475	if((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1)))
1476	& then
1477	iref=isoil
1478	exit
1479	endif
1480	enddo
1481	! endif
1482
1483	! write(,)'iref:',iref,' jref:',jref
1484	! write(,)'layer',layer
1485	! write(,)'mlayer',mlayer
1486
1487	! thermal inertia of the ice:
1488	ierr=1
1489	do while (ierr.ne.0)
1490	write(,)'What is the value of subterranean ice thermal inert
1491	&ia? (e.g.: 2000)'
1492	read(,,iostat=ierr)iceith
1493	enddo ! of do while
1494
1495	! recast ithfi() onto ith()
1496	call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
1497
1498	do j=1,jref
1499	! write(,)'j:',j,'rlatu(j)180./pi:',rlatu(j)180./pi
1500	do i=1,iip1 ! loop on longitudes
1501	! Build "equivalent" thermal inertia for the mixed layer
1502	ith(i,j,iref+1)=sqrt((layer(iref+1)-layer(iref))/
1503	& (((val2-layer(iref))/(ith(i,j,iref)**2))+
1504	& ((layer(iref+1)-val2)/(iceith)**2)))
1505	! Set thermal inertia of lower layers
1506	do isoil=iref+2,nsoilmx
1507	ith(i,j,isoil)=iceith ! ice
1508	enddo
1509	enddo ! of do i=1,iip1
1510	enddo ! of do j=1,jjp1
1511
1512
1513	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
1514
1515	! do i=1,nsoilmx
1516	! write(,)'i:',i,'ithfi(1,i):',ithfi(1,i)
1517	! enddo
1518
1519
1520	! subsoilice_s: Put deep ice layer in southern hemisphere soil
1521	! ------------------------------------------------------------
1522
1523	else if (trim(modif).eq.'subsoilice_s') then
1524
1525	write(,)'From which latitude (in deg.), down to the south pol
1526	&e, should we put subterranean ice?'
1527	ierr=1
1528	do while (ierr.ne.0)
1529	read(,,iostat=ierr) val
1530	if (ierr.eq.0) then ! got a value
1531	! do a sanity check
1532	if((val.gt.0.).or.(val.lt.-90)) then
1533	write(,)'Latitude should be between 0 and -90 deg. !!!'
1534	ierr=1
1535	else ! find corresponding jref (nearest latitude)
1536	! note: rlatu(:) contains decreasing values of latitude
1537	! starting from PI/2 to -PI/2
1538	do j=1,jjp1
1539	if ((rlatu(j)*180./pi.ge.val).and.
1540	& (rlatu(j+1)*180./pi.le.val)) then
1541	! find which grid point is nearest to val:
1542	if (abs(rlatu(j)*180./pi-val).le.
1543	& abs((rlatu(j+1)*180./pi-val))) then
1544	jref=j
1545	else
1546	jref=j+1
1547	endif
1548
1549	write(,)'Will use nearest grid latitude which is:',
1550	& rlatu(jref)*180./pi
1551	endif
1552	enddo ! of do j=1,jjp1
1553	endif ! of if((val.lt.0.).or.(val.gt.90))
1554	endif !of if (ierr.eq.0)
1555	enddo ! of do while
1556
1557	! Build layers() (as in soil_settings.F)
1558	val2=sqrt(mlayer(0)*mlayer(1))
1559	val3=mlayer(1)/mlayer(0)
1560	do isoil=1,nsoilmx
1561	layer(isoil)=val2(val3*(isoil-1))
1562	enddo
1563
1564	write(,)'At which depth (in m.) does the ice layer begin?'
1565	write(,)'(currently, the deepest soil layer extends down to:'
1566	& ,layer(nsoilmx),')'
1567	ierr=1
1568	do while (ierr.ne.0)
1569	read(,,iostat=ierr) val2
1570	! write(,)'val2:',val2,'ierr=',ierr
1571	if (ierr.eq.0) then ! got a value, but do a sanity check
1572	if(val2.gt.layer(nsoilmx)) then
1573	write(,)'Depth should be less than ',layer(nsoilmx)
1574	ierr=1
1575	endif
1576	if(val2.lt.layer(1)) then
1577	write(,)'Depth should be more than ',layer(1)
1578	ierr=1
1579	endif
1580	endif
1581	enddo ! of do while
1582
1583	! find the reference index iref the depth corresponds to
1584	do isoil=1,nsoilmx-1
1585	if((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1)))
1586	& then
1587	iref=isoil
1588	exit
1589	endif
1590	enddo
1591
1592	! write(,)'iref:',iref,' jref:',jref
1593
1594	! thermal inertia of the ice:
1595	ierr=1
1596	do while (ierr.ne.0)
1597	write(,)'What is the value of subterranean ice thermal inert
1598	&ia? (e.g.: 2000)'
1599	read(,,iostat=ierr)iceith
1600	enddo ! of do while
1601
1602	! recast ithfi() onto ith()
1603	call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
1604
1605	do j=jref,jjp1
1606	! write(,)'j:',j,'rlatu(j)180./pi:',rlatu(j)180./pi
1607	do i=1,iip1 ! loop on longitudes
1608	! Build "equivalent" thermal inertia for the mixed layer
1609	ith(i,j,iref+1)=sqrt((layer(iref+1)-layer(iref))/
1610	& (((val2-layer(iref))/(ith(i,j,iref)**2))+
1611	& ((layer(iref+1)-val2)/(iceith)**2)))
1612	! Set thermal inertia of lower layers
1613	do isoil=iref+2,nsoilmx
1614	ith(i,j,isoil)=iceith ! ice
1615	enddo
1616	enddo ! of do i=1,iip1
1617	enddo ! of do j=jref,jjp1
1618
1619
1620	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
1621
1622	c 'mons_ice' : use MONS data to build subsurface ice table
1623	c --------------------------------------------------------
1624	else if (trim(modif).eq.'mons_ice') then
1625
1626	! 1. Load MONS data
1627	call load_MONS_data(MONS_Hdn,MONS_d21)
1628
1629	! 2. Get parameters from user
1630	ierr=1
1631	do while (ierr.ne.0)
1632	write(,) "Coefficient to apply to MONS 'depth' in Northern",
1633	& " Hemisphere?"
1634	write(,) " (should be somewhere between 3.2e-4 and 1.3e-3)"
1635	read(,,iostat=ierr) MONS_coeffN
1636	enddo
1637	ierr=1
1638	do while (ierr.ne.0)
1639	write(,) "Coefficient to apply to MONS 'depth' in Southern",
1640	& " Hemisphere?"
1641	write(,) " (should be somewhere between 3.2e-4 and 1.3e-3)"
1642	read(,,iostat=ierr) MONS_coeffS
1643	enddo
1644	ierr=1
1645	do while (ierr.ne.0)
1646	write(,) "Value of subterranean ice thermal inertia ",
1647	& " in Northern hemisphere?"
1648	write(,) " (e.g.: 2000, or perhaps 2290)"
1649	! read(,,iostat=ierr) iceith
1650	read(,,iostat=ierr) iceithN
1651	enddo
1652	ierr=1
1653	do while (ierr.ne.0)
1654	write(,) "Value of subterranean ice thermal inertia ",
1655	& " in Southern hemisphere?"
1656	write(,) " (e.g.: 2000, or perhaps 2290)"
1657	! read(,,iostat=ierr) iceith
1658	read(,,iostat=ierr) iceithS
1659	enddo
1660
1661	! 3. Build subterranean thermal inertia
1662
1663	! initialise subsurface inertia with reference surface values
1664	do isoil=1,nsoilmx
1665	ithfi(1:ngridmx,isoil)=surfithfi(1:ngridmx)
1666	enddo
1667	! recast ithfi() onto ith()
1668	call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
1669
1670	do i=1,iip1 ! loop on longitudes
1671	do j=1,jjp1 ! loop on latitudes
1672	! set MONS_coeff
1673	if (rlatu(j).ge.0) then ! northern hemisphere
1674	! N.B: rlatu(:) contains decreasing values of latitude
1675	! starting from PI/2 to -PI/2
1676	MONS_coeff=MONS_coeffN
1677	iceith=iceithN
1678	else ! southern hemisphere
1679	MONS_coeff=MONS_coeffS
1680	iceith=iceithS
1681	endif
1682	! check if we should put subterranean ice
1683	if (MONS_Hdn(i,j).ge.14.0) then ! no ice if Hdn<14%
1684	! compute depth at which ice lies:
1685	val=MONS_d21(i,j)*MONS_coeff
1686	! compute val2= the diurnal skin depth of surface inertia
1687	! assuming a volumetric heat cap. of C=1.e6 J.m-3.K-1
1688	val2=ith(i,j,1)1.e-6sqrt(88775./3.14159)
1689	if (val.lt.val2) then
1690	! ice must be below the (surface inertia) diurnal skin depth
1691	val=val2
1692	endif
1693	if (val.lt.layer(nsoilmx)) then ! subterranean ice
1694	! find the reference index iref that depth corresponds to
1695	iref=0
1696	do isoil=1,nsoilmx-1
1697	if ((val.ge.layer(isoil)).and.(val.lt.layer(isoil+1)))
1698	& then
1699	iref=isoil
1700	exit
1701	endif
1702	enddo
1703	! Build "equivalent" thermal inertia for the mixed layer
1704	ith(i,j,iref+1)=sqrt((layer(iref+1)-layer(iref))/
1705	& (((val-layer(iref))/(ith(i,j,iref+1)**2))+
1706	& ((layer(iref+1)-val)/(iceith)**2)))
1707	! Set thermal inertia of lower layers
1708	do isoil=iref+2,nsoilmx
1709	ith(i,j,isoil)=iceith
1710	enddo
1711	endif ! of if (val.lt.layer(nsoilmx))
1712	endif ! of if (MONS_Hdn(i,j).lt.14.0)
1713	enddo ! do j=1,jjp1
1714	enddo ! do i=1,iip1
1715
1716	! Check:
1717	! do i=1,iip1
1718	! do j=1,jjp1
1719	! do isoil=1,nsoilmx
1720	! write(77,*) i,j,isoil," ",ith(i,j,isoil)
1721	! enddo
1722	! enddo
1723	! enddo
1724
1725	! recast ith() into ithfi()
1726	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
1727
1728	else if (trim(modif) .eq. 'nslope') then
1729	write(,) 'set a new number of subgrid scale slope'
1730	write(,) 'Current value=', nslope
1731	write(,) 'Enter value for nslope (ex: 1,5,7)?'
1732	ierr=1
1733	do while (ierr.ne.0)
1734	read(,,iostat=ierr) nslope_new
1735	enddo
1736
1737	write(,) 'This can take some time...'
1738	write(,) 'You can go grab a coffee and relax a bit'
1739
1740	if(nslope.eq.nslope_new) then
1741	write(,) 'The number of subslope you entered is the same'
1742	write(,) 'as the number written in startfi.nc. '
1743	write(,) 'Nothing will be done'
1744	else
1745
1746	nslope_old=nslope
1747	nslope=nslope_new
1748
1749	call end_comslope_h
1750	call ini_comslope_h(ngridmx,nslope_new)
1751
1752	allocate(default_def_slope(nslope_new+1))
1753	!Sub-grid scale subslopes
1754	if (nslope_new.eq.7) then
1755	default_def_slope(1) = -43.
1756	default_def_slope(2) = -19.
1757	default_def_slope(3) = -9.
1758	default_def_slope(4) = -3.
1759	default_def_slope(5) = 3.
1760	default_def_slope(6) = 9.
1761	default_def_slope(7) = 19.
1762	default_def_slope(8) = 43.
1763	elseif (nslope_new.eq.5) then
1764	default_def_slope(1) = -43.
1765	default_def_slope(2) = -9.
1766	default_def_slope(3) = -3.
1767	default_def_slope(4) = 3.
1768	default_def_slope(5) = 9.
1769	default_def_slope(6) = 43.
1770	elseif (nslope_new.eq.1) then
1771	default_def_slope(1) = -50.
1772	default_def_slope(2) = 50.
1773	endif
1774
1775	do islope=1,nslope_new+1
1776	def_slope(islope) = default_def_slope(islope)
1777	enddo
1778
1779	do islope=1,nslope_new
1780	def_slope_mean(islope) =(def_slope(islope)+def_slope(islope+1))/2.
1781	enddo
1782
1783	iflat = 1
1784	DO islope=2,nslope_new
1785	IF(abs(def_slope_mean(islope)).lt.
1786	& abs(def_slope_mean(iflat)))THEN
1787	iflat = islope
1788	ENDIF
1789	ENDDO
1790
1791	call subslope_mola(ngridmx,nslope_new,def_slope,subslope_dist)
1792
1793	! Surfdat related stuff
1794
1795	allocate(tsurf_old_slope(ngridmx,nslope_old))
1796	allocate(qsurf_old_slope(ngridmx,nqtot,nslope_old))
1797	allocate(emis_old_slope(ngridmx,nslope_old))
1798	allocate(watercap_old_slope(ngridmx,nslope_old))
1799
1800	tsurf_old_slope(:,:)=tsurf(:,:)
1801	qsurf_old_slope(:,:,:)=qsurf(:,:,:)
1802	emis_old_slope(:,:)=emis(:,:)
1803	watercap_old_slope(:,:)=watercap(:,:)
1804
1805	call end_surfdat_h_slope_var
1806	call ini_surfdat_h_slope_var(ngridmx,nqtot,nslope_new)
1807
1808	! Comsoil related stuff (tsoil)
1809
1810	allocate(tsoil_old_slope(ngridmx,nsoilmx,nslope_old))
1811
1812	tsoil_old_slope(:,:,:)=tsoil(:,:,:)
1813
1814	call end_comsoil_h_slope_var
1815	call ini_comsoil_h_slope_var(ngridmx,nslope_new)
1816
1817	! Dimradmars related stuff (albedo)
1818
1819	allocate(albedo_old_slope(ngridmx,2,nslope_old))
1820	albedo_old_slope(:,:,:)=albedo(:,:,:)
1821
1822	call end_dimradmars_mod_slope_var
1823	call ini_dimradmars_mod_slope_var(ngridmx,nslope_new)
1824
1825	if(nslope_old.eq.1 .and. nslope_new.gt.1) then
1826	do islope=1,nslope_new
1827	tsurf(:,islope)=tsurf_old_slope(:,1)
1828	qsurf(:,:,islope)=qsurf_old_slope(:,:,1)
1829	emis(:,islope)=emis_old_slope(:,1)
1830	watercap(:,islope)=watercap_old_slope(:,1)
1831	tsoil(:,:,islope)=tsoil_old_slope(:,:,1)
1832	albedo(:,:,islope)=albedo_old_slope(:,:,1)
1833	enddo
1834	elseif(nslope_new.eq.1) then
1835	tsurf(:,1)=tsurf_old_slope(:,iflat)
1836	qsurf(:,:,1)=qsurf_old_slope(:,:,iflat)
1837	emis(:,1)=emis_old_slope(:,iflat)
1838	watercap(:,1)=watercap_old_slope(:,iflat)
1839	tsoil(:,:,1)=tsoil_old_slope(:,:,iflat)
1840	albedo(:,:,1)=albedo_old_slope(:,:,iflat)
1841	elseif(nslope_old.eq.5 .and. nslope_new.eq.7) then
1842	do islope=1,nslope_new
1843	tsurf(:,islope)=tsurf_old_slope(:,iflat)
1844	qsurf(:,:,islope)=qsurf_old_slope(:,:,iflat)
1845	emis(:,islope)=emis_old_slope(:,iflat)
1846	watercap(:,islope)=watercap_old_slope(:,iflat)
1847	tsoil(:,:,islope)=tsoil_old_slope(:,:,iflat)
1848	albedo(:,:,islope)=albedo_old_slope(:,:,iflat)
1849	enddo
1850	elseif(nslope_old.eq.7 .and. nslope_new.eq.5) then
1851	do islope=1,nslope_new
1852	tsurf(:,islope)=tsurf_old_slope(:,iflat)
1853	qsurf(:,:,islope)=qsurf_old_slope(:,:,iflat)
1854	emis(:,islope)=emis_old_slope(:,iflat)
1855	watercap(:,islope)=watercap_old_slope(:,iflat)
1856	tsoil(:,:,islope)=tsoil_old_slope(:,:,iflat)
1857	albedo(:,:,islope)=albedo_old_slope(:,:,iflat)
1858	enddo
1859	else
1860	write(,)' Problem choice of nslope'
1861	write(,)' Value not taken into account'
1862	CALL ABORT
1863	endif
1864
1865	endif !nslope=nslope_new
1866
1867	else
1868	write(,) ' Unknown (misspelled?) option!!!'
1869	end if ! of if (trim(modif) .eq. '...') elseif ...
1870
1871	enddo ! of do ! infinite loop on liste of changes
1872
1873	999 continue
1874
1875
1876	c=======================================================================
1877	c Correct pressure on the new grid (menu 0)
1878	c=======================================================================
1879
1880	if (choix_1.eq.0) then
1881	r = 1000.*8.31/mugaz
1882
1883	do j=1,jjp1
1884	do i=1,iip1
1885	ps(i,j) = ps(i,j) *
1886	. exp((phisold_newgrid(i,j)-phis(i,j)) /
1887	. (t(i,j,1) * r))
1888	end do
1889	end do
1890
1891	c periodicity of surface ps in longitude
1892	do j=1,jjp1
1893	ps(1,j) = ps(iip1,j)
1894	end do
1895	end if
1896
1897	c=======================================================================
1898	c=======================================================================
1899
1900	c=======================================================================
1901	c Initialisation de la physique / ecriture de newstartfi :
1902	c=======================================================================
1903
1904
1905	CALL inifilr
1906	CALL pression(ip1jmp1, ap, bp, ps, p3d)
1907
1908	c-----------------------------------------------------------------------
1909	c Initialisation pks:
1910	c-----------------------------------------------------------------------
1911
1912	CALL exner_hyb(ip1jmp1, ps, p3d, pks, pk, pkf)
1913	! Calcul de la temperature potentielle teta
1914
1915	if (flagiso) then
1916	DO l=1,llm
1917	DO j=1,jjp1
1918	DO i=1,iim
1919	teta(i,j,l) = Tiso * cpp/pk(i,j,l)
1920	ENDDO
1921	teta (iip1,j,l)= teta (1,j,l)
1922	ENDDO
1923	ENDDO
1924	else if (choix_1.eq.0) then
1925	DO l=1,llm
1926	DO j=1,jjp1
1927	DO i=1,iim
1928	teta(i,j,l) = t(i,j,l) * cpp/pk(i,j,l)
1929	ENDDO
1930	teta (iip1,j,l)= teta (1,j,l)
1931	ENDDO
1932	ENDDO
1933	endif
1934
1935	C Calcul intermediaire
1936	c
1937	if (choix_1.eq.0) then
1938	CALL massdair( p3d, masse )
1939	c
1940	print *,' ALPHAX ',alphax
1941
1942	DO l = 1, llm
1943	DO i = 1, iim
1944	xppn(i) = aire( i, 1 ) * masse( i , 1 , l )
1945	xpps(i) = aire( i,jjp1 ) * masse( i , jjp1 , l )
1946	ENDDO
1947	xpn = SUM(xppn)/apoln
1948	xps = SUM(xpps)/apols
1949	DO i = 1, iip1
1950	masse( i , 1 , l ) = xpn
1951	masse( i , jjp1 , l ) = xps
1952	ENDDO
1953	ENDDO
1954	endif
1955	phis(iip1,:) = phis(1,:)
1956
1957	c CALL inidissip ( lstardis, nitergdiv, nitergrot, niterh,
1958	c * tetagdiv, tetagrot , tetatemp )
1959	itau=0
1960	if (choix_1.eq.0) then
1961	day_ini=int(date)
1962	hour_ini=date-int(date)
1963	endif
1964	c
1965	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
1966
1967	CALL caldyn0( itau,ucov,vcov,teta,ps,masse,pk,phis ,
1968	* phi,w, pbaru,pbarv,day_ini+time )
1969	c CALL caldyn
1970	c $ ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
1971	c $ phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, day_ini )
1972
1973	CALL dynredem0("restart.nc",day_ini,phis)
1974	CALL dynredem1("restart.nc",hour_ini,vcov,ucov,teta,q,
1975	. masse,ps)
1976	C
1977	C Ecriture etat initial physique
1978	C
1979	call physdem0("restartfi.nc",longitude,latitude,
1980	& nsoilmx,ngridmx,llm,
1981	& nqtot,dtphys,real(day_ini),0.0,cell_area,
1982	& albfi,ithfi,def_slope,
1983	& subslope_dist)
1984	call physdem1("restartfi.nc",nsoilmx,ngridmx,llm,nqtot,
1985	& dtphys,hour_ini,
1986	& tsurf,tsoil,albedo,emis,q2,qsurf,tauscaling,
1987	& totcloudfrac,wstar,watercap)
1988
1989	c=======================================================================
1990	c Formats
1991	c=======================================================================
1992
1993	1 FORMAT(//10x,'la valeur de im =',i4,2x,'lue sur le fichier de dema
1994	*rrage est differente de la valeur parametree iim =',i4//)
1995	2 FORMAT(//10x,'la valeur de jm =',i4,2x,'lue sur le fichier de dema
1996	*rrage est differente de la valeur parametree jjm =',i4//)
1997	3 FORMAT(//10x,'la valeur de lllm =',i4,2x,'lue sur le fichier demar
1998	*rage est differente de la valeur parametree llm =',i4//)
1999
2000	write(,) "newstart: All is well that ends well."
2001
2002	end
2003
2004	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2005	subroutine load_MONS_data(MONS_Hdn,MONS_d21)
2006
2007	use datafile_mod, only:datadir
2008
2009	implicit none
2010	! routine to load Benedicte Diez MONS dataset, fill in date in southern
2011	! polar region, and interpolate the result onto the GCM grid
2012	include"dimensions.h"
2013	include"paramet.h"
2014	include"comgeom.h"
2015	! arguments:
2016	real,intent(out) :: MONS_Hdn(iip1,jjp1) ! Hdn: %WEH=Mass fraction of H2O
2017	real,intent(out) :: MONS_d21(iip1,jjp1) ! ice table "depth" (in kg/m2)
2018	! N.B MONS datasets should be of dimension (iip1,jjp1)
2019	! local variables:
2020	character(len=88) :: filename="results_MONS_lat_lon_H_depth.txt"
2021	character(len=88) :: txt ! to store some text
2022	integer :: ierr,i,j
2023	integer,parameter :: nblon=180 ! number of longitudes of MONS datasets
2024	integer,parameter :: nblat=90 ! number of latitudes of MONS datasets
2025	real :: pi
2026	real :: longitudes(nblon) ! MONS dataset longitudes
2027	real :: latitudes(nblat) ! MONS dataset latitudes
2028	! MONS dataset: mass fraction of H2O where H is assumed to be in H2O
2029	real :: Hdn(nblon,nblat)
2030	real :: d21(nblon,nblat)! MONS dataset "depth" (g/cm2)
2031
2032	! Extended MONS dataset (for interp_horiz)
2033	real :: Hdnx(nblon+1,nblat)
2034	real :: d21x(nblon+1,nblat)
2035	real :: lon_bound(nblon+1) ! longitude boundaries
2036	real :: lat_bound(nblat-1) ! latitude boundaries
2037
2038	! 1. Initializations:
2039
2040	write(,) "Loading MONS data"
2041
2042	! Open MONS datafile:
2043	open(42,file=trim(datadir)//"/"//trim(filename),
2044	& status="old",iostat=ierr)
2045	if (ierr/=0) then
2046	write(,) "Error in load_MONS_data:"
2047	write(,) "Failed opening file ",
2048	& trim(datadir)//"/"//trim(filename)
2049	write(,)'1) You can change this directory address in ',
2050	& 'callfis.def with'
2051	write(,)' datadir=/path/to/datafiles'
2052	write(,)'2) If necessary ',trim(filename),
2053	& ' (and other datafiles)'
2054	write(,)' can be obtained online at:'
2055	write(,)'http://www.lmd.jussieu.fr/~lmdz/planets/mars/datadir'
2056	CALL ABORT
2057	else ! skip first line of file (dummy read)
2058	read(42,*) txt
2059	endif
2060
2061	pi=2.*asin(1.)
2062
2063	!2. Load MONS data (on MONS grid)
2064	do j=1,nblat
2065	do i=1,nblon
2066	! swap latitude index so latitudes go from north pole to south pole:
2067	read(42,*) latitudes(nblat-j+1),longitudes(i),
2068	& Hdn(i,nblat-j+1),d21(i,nblat-j+1)
2069	! multiply d21 by 10 to convert from g/cm2 to kg/m2
2070	d21(i,nblat-j+1)=d21(i,nblat-j+1)*10.0
2071	enddo
2072	enddo
2073	close(42)
2074
2075	! there is unfortunately no d21 data for latitudes -77 to -90
2076	! so we build some by linear interpolation between values at -75
2077	! and assuming d21=0 at the pole
2078	do j=84,90 ! latitudes(84)=-77 ; latitudes(83)=-75
2079	do i=1,nblon
2080	d21(i,j)=d21(i,83)*((latitudes(j)+90)/15.0)
2081	enddo
2082	enddo
2083
2084	! 3. Build extended MONS dataset & boundaries (for interp_horiz)
2085	! longitude boundaries (in radians):
2086	do i=1,nblon
2087	! NB: MONS data is every 2 degrees in longitude
2088	lon_bound(i)=(longitudes(i)+1.0)*pi/180.0
2089	enddo
2090	! extra 'modulo' value
2091	lon_bound(nblon+1)=lon_bound(1)+2.0*pi
2092
2093	! latitude boundaries (in radians):
2094	do j=1,nblat-1
2095	! NB: Mons data is every 2 degrees in latitude
2096	lat_bound(j)=(latitudes(j)-1.0)*pi/180.0
2097	enddo
2098
2099	! MONS datasets:
2100	do j=1,nblat
2101	Hdnx(1:nblon,j)=Hdn(1:nblon,j)
2102	Hdnx(nblon+1,j)=Hdnx(1,j)
2103	d21x(1:nblon,j)=d21(1:nblon,j)
2104	d21x(nblon+1,j)=d21x(1,j)
2105	enddo
2106
2107	! Interpolate onto GCM grid
2108	call interp_horiz(Hdnx,MONS_Hdn,nblon,nblat-1,iim,jjm,1,
2109	& lon_bound,lat_bound,rlonu,rlatv)
2110	call interp_horiz(d21x,MONS_d21,nblon,nblat-1,iim,jjm,1,
2111	& lon_bound,lat_bound,rlonu,rlatv)
2112
2113	end subroutine

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

Download in other formats:

Original Format