Context Navigation

newstart.F @ 433

Last change on this file since 433 was 321, checked in by emillour, 14 years ago

Mars GCM:

Corrected small bug in newstart: initracer was not always used and thus some tracer indexes (igm_co2, igcm_h2o_vap,...) were not set. This however means that we now also call inifis from newstart and that we read in flags set in 'callphys.def' (required for sanity checks in initracer). Also adapted 'inichim_readcallphys': removed some obsolescent tests on number of tracers for given combinations of options.

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

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

Download in other formats:

Original Format