Context Navigation

newstart.F @ 2551

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

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

Download in other formats:

Original Format