Context Navigation

newstart.F @ 676

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

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

Download in other formats:

Original Format