Context Navigation

newstart.F @ 1104

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

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

Download in other formats:

Original Format