Context Navigation

newstart.F @ 3316

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

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

Download in other formats:

Original Format