Context Navigation

newstart.F @ 563

Last change on this file since 563 was 563, checked in by emillour, 13 years ago

Mars GCM:

Added option "q=profile" in newstart to initialize a given tracer with

a profile specified in an ASCII file.

Minor fix on lect_start_archive : when tracers were not found, a default

value was asked to user, but twice! Once is enough.

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

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

Download in other formats:

Original Format