Context Navigation

newstart.F @ 1130

Last change on this file since 1130 was 1130, checked in by emillour, 11 years ago

Mars GCM:
Series of changes to enable running in parallel (using LMDZ.COMMON dynamics);
Current LMDZ.MARS can still notheless be compiled and run in serial mode
"as previously".
Summary of main changes:

Main programs (newstart, start2archive, xvik) that used to be in dyn3d have been moved to phymars.
dyn3d/control.h is now module control_mod.F90
rearanged input/outputs routines everywhere to handle serial/MPI cases. physdem.F => phyredem.F90 , phyetat0.F => phyetat0.F90 ; all read/write routines for startfi files are gathered in module iostart.F90
added parallelism related routines init_phys_lmdz.F90, comgeomphy.F90, dimphy.F90, iniphysiq.F90, mod_grid_phy_lmdz.F90, mod_phys_lmdz_mpi_data.F90, mod_phys_lmdz_mpi_transfert.F90, mod_phys_lmdz_omp_data.F90, mod_phys_lmdz_omp_transfert.F90, mod_phys_lmdz_para.F90, mod_phys_lmdz_transfert_para.F90 in phymars and mod_const_mpi.F90 in dyn3d (for compliance with parallel case)
created generic routines 'planetwide_maxval' and 'planetwide_minval', in module "planetwide_mod", that enable obtaining the min and max of a field over the whole planet.

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

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

Download in other formats:

Original Format