Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

newstart.F.tbold @ 3436

Last change on this file since 3436 was 3175, checked in by emillour, 11 months ago
Pluto PCM: Add the old Pluto LMDZ for reference (required prior step to making an LMDZ.PLUTO using the same framework as the other physics packages). TB+EM
File size: 66.0 KB

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

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

Download in other formats:

Original Format