Context Navigation

newstart.F @ 205

Last change on this file since 205 was 135, checked in by aslmd, 14 years ago
CHANGEMENT ARBORESCENCE ETAPE 2 -- NON COMPLET
File size: 53.8 KB

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

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

Download in other formats:

Original Format