Context Navigation

newstart.F @ 3302

Last change on this file since 3302 was 3203, checked in by jbclement, 11 months ago

Mars PCM:

Addition of the possibility to use subslopes parametization in 1D. To do so, put 'nslope=x' in "run.def" where x (1, 3, 5 or 7) is the number of slopes you want to, or modify it in an appropriate "startfi.nc". Then, default subslopes definition and distribution are used by the model. The already existing case of using one slope whose inclination and orientation are defined in "run.def" is still possible.
Some cleanings throughout the code, in particular related to unused variables.

JBC

File size: 76.8 KB

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

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

Download in other formats:

Original Format