Context Navigation

newstart.F @ 300

Last change on this file since 300 was 253, checked in by emillour, 14 years ago

Generic GCM

Massive update to version 0.7

EM+RW

File size: 65.4 KB

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

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

Download in other formats:

Original Format