Context Navigation

newstart.F @ 1127

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

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

Download in other formats:

Original Format