Context Navigation

newstart.F @ 1212

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

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

Download in other formats:

Original Format