Context Navigation

newstart.F @ 3604

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

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

Download in other formats:

Original Format