Context Navigation

newstart.F @ 699

Last change on this file since 699 was 699, checked in by emillour, 13 years ago

Generic GCM:

Corrected the polar mesh surface area which was wrong in the physics (changes in phyetat0.F, calfis.F and newstart.F)
Some cleanup in newstart.F (removed some obsolete "Mars" options: mons_ice,..) and also added option "q=profile" to initialize a tracer with a profile read from file "profile_tracername"

File size: 47.9 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	real :: profile(llm+1) ! to store an atmospheric profile + surface value
148
149	! added by RW for test
150	real pmean, phi0
151
152	! added by BC for equilibrium temperature startup
153	real teque
154
155	! added by BC for cloud fraction setup
156	REAL hice(ngridmx),cloudfrac(ngridmx,nlayermx)
157	REAL totalfrac(ngridmx)
158
159	! added by RW for nuketharsis
160	real fact1
161	real fact2
162
163
164	c sortie visu pour les champs dynamiques
165	c---------------------------------------
166	! INTEGER :: visuid
167	! real :: time_step,t_ops,t_wrt
168	! CHARACTER*80 :: visu_file
169
170	cpp = 0.
171	preff = 0.
172	pa = 0. ! to ensure disaster rather than minor error if we don`t
173	! make deliberate choice of these values elsewhere.
174
175	c=======================================================================
176	c Choice of the start file(s) to use
177	c=======================================================================
178
179	write(,) 'From which kind of files do you want to create new',
180	. 'start and startfi files'
181	write(,) ' 0 - from a file start_archive'
182	write(,) ' 1 - from files start and startfi'
183
184	c-----------------------------------------------------------------------
185	c Open file(s) to modify (start or start_archive)
186	c-----------------------------------------------------------------------
187
188	DO
189	read(,,iostat=ierr) choix_1
190	if ((choix_1 /= 0).OR.(choix_1 /=1)) EXIT
191	ENDDO
192
193	c Open start_archive
194	c ~~~~~~~~~~~~~~~~~~~~~~~~~~
195	if (choix_1.eq.0) then
196
197	write(,) 'Creating start files from:'
198	write(,) './start_archive.nc'
199	write(,)
200	fichnom = 'start_archive.nc'
201	ierr = NF_OPEN (fichnom, NF_NOWRITE,nid)
202	IF (ierr.NE.NF_NOERR) THEN
203	write(6,*)' Problem opening file:',fichnom
204	write(6,*)' ierr = ', ierr
205	CALL ABORT
206	ENDIF
207	tab0 = 50
208	Lmodif = 1
209
210	c OR open start and startfi files
211	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
212	else
213	write(,) 'Creating start files from:'
214	write(,) './start.nc and ./startfi.nc'
215	write(,)
216	fichnom = 'start.nc'
217	ierr = NF_OPEN (fichnom, NF_NOWRITE,nid_dyn)
218	IF (ierr.NE.NF_NOERR) THEN
219	write(6,*)' Problem opening file:',fichnom
220	write(6,*)' ierr = ', ierr
221	CALL ABORT
222	ENDIF
223
224	fichnom = 'startfi.nc'
225	ierr = NF_OPEN (fichnom, NF_NOWRITE,nid_fi)
226	IF (ierr.NE.NF_NOERR) THEN
227	write(6,*)' Problem opening file:',fichnom
228	write(6,*)' ierr = ', ierr
229	CALL ABORT
230	ENDIF
231
232	tab0 = 0
233	Lmodif = 0
234
235	endif
236
237
238	c=======================================================================
239	c INITIALISATIONS DIVERSES
240	c=======================================================================
241	! Load tracer names:
242	call iniadvtrac(nq,numvanle)
243	! tnom(:) now contains tracer names
244	! JL12 we will need the tracer names to read start in dyneta0
245
246	c-----------------------------------------------------------------------
247	c Lecture du tableau des parametres du run (pour la dynamique)
248	c-----------------------------------------------------------------------
249
250	if (choix_1.eq.0) then
251
252	write(,) 'reading tab_cntrl START_ARCHIVE'
253	c
254	ierr = NF_INQ_VARID (nid, "controle", nvarid)
255	#ifdef NC_DOUBLE
256	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tab_cntrl)
257	#else
258	ierr = NF_GET_VAR_REAL(nid, nvarid, tab_cntrl)
259	#endif
260	c
261	else if (choix_1.eq.1) then
262
263	write(,) 'reading tab_cntrl START'
264	c
265	ierr = NF_INQ_VARID (nid_dyn, "controle", nvarid)
266	#ifdef NC_DOUBLE
267	ierr = NF_GET_VAR_DOUBLE(nid_dyn, nvarid, tab_cntrl)
268	#else
269	ierr = NF_GET_VAR_REAL(nid_dyn, nvarid, tab_cntrl)
270	#endif
271	c
272	write(,) 'reading tab_cntrl STARTFI'
273	c
274	ierr = NF_INQ_VARID (nid_fi, "controle", nvarid)
275	#ifdef NC_DOUBLE
276	ierr = NF_GET_VAR_DOUBLE(nid_fi, nvarid, tab_cntrl_bis)
277	#else
278	ierr = NF_GET_VAR_REAL(nid_fi, nvarid, tab_cntrl_bis)
279	#endif
280	c
281	do i=1,50
282	tab_cntrl(i+50)=tab_cntrl_bis(i)
283	enddo
284	write(,) 'printing tab_cntrl', tab_cntrl
285	do i=1,100
286	write(,) i,tab_cntrl(i)
287	enddo
288
289	! Lmodif set to 0 to disable modifications possibility in phyeta0
290	write(,) 'Reading file START'
291	fichnom = 'start.nc'
292	CALL dynetat0(fichnom,nqmx,vcov,ucov,teta,q,masse,
293	. ps,phis,time)
294
295	write(,) 'Reading file STARTFI'
296	fichnom = 'startfi.nc'
297	CALL phyetat0 (fichnom,tab0,Lmodif,nsoilmx,nqmx,
298	. day_ini,time,
299	. tsurf,tsoil,emis,q2,qsurf, !) ! temporary modif by RDW
300	. cloudfrac,totalfrac,hice)
301
302	! copy albedo and soil thermal inertia
303	do i=1,ngridmx
304	albfi(i) = albedodat(i)
305	do j=1,nsoilmx
306	ithfi(i,j) = inertiedat(i,j)
307	enddo
308	! build a surfithfi(:) using 1st layer of ithfi(:), which might
309	! be neede later on if reinitializing soil thermal inertia
310	surfithfi(i)=ithfi(i,1)
311	enddo
312
313
314	endif
315	c-----------------------------------------------------------------------
316	c Initialisation des constantes dynamique
317	c-----------------------------------------------------------------------
318
319	kappa = tab_cntrl(9)
320	etot0 = tab_cntrl(12)
321	ptot0 = tab_cntrl(13)
322	ztot0 = tab_cntrl(14)
323	stot0 = tab_cntrl(15)
324	ang0 = tab_cntrl(16)
325	write(,) "Newstart: kappa,etot0,ptot0,ztot0,stot0,ang0"
326	write(,) kappa,etot0,ptot0,ztot0,stot0,ang0
327
328	! for vertical coordinate
329	preff=tab_cntrl(18) ! reference surface pressure
330	pa=tab_cntrl(17) ! reference pressure at which coord is purely pressure
331	!NB: in start_archive files tab_cntrl(17)=tab_cntrl(18)=0
332	write(,) "Newstart: preff=",preff," pa=",pa
333	yes=' '
334	do while ((yes.ne.'y').and.(yes.ne.'n'))
335	write(,) "Change the values of preff and pa? (y/n)"
336	read(*,fmt='(a)') yes
337	end do
338	if (yes.eq.'y') then
339	write(,)"New value of reference surface pressure preff?"
340	write(,)" (for Mars, typically preff=610)"
341	read(,) preff
342	write(,)"New value of reference pressure pa for purely"
343	write(,)"pressure levels (for hybrid coordinates)?"
344	write(,)" (for Mars, typically pa=20)"
345	read(,) pa
346	endif
347	c-----------------------------------------------------------------------
348	c Lecture du tab_cntrl et initialisation des constantes physiques
349	c - pour start: Lmodif = 0 => pas de modifications possibles
350	c (modif dans le tabfi de readfi + loin)
351	c - pour start_archive: Lmodif = 1 => modifications possibles
352	c-----------------------------------------------------------------------
353	if (choix_1.eq.0) then
354	call tabfi (nid,Lmodif,tab0,day_ini,lllm,p_rad,
355	. p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
356	else if (choix_1.eq.1) then
357	Lmodif=1 ! Lmodif set to 1 to allow modifications in phyeta0
358	call tabfi (nid_fi,Lmodif,tab0,day_ini,lllm,p_rad,
359	. p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
360	endif
361
362	rad = p_rad
363	omeg = p_omeg
364	g = p_g
365	cpp = p_cpp
366	mugaz = p_mugaz
367	daysec = p_daysec
368
369	kappa = 8.3141000./(p_mugazp_cpp) ! added by RDW
370
371
372	c=======================================================================
373	c INITIALISATIONS DIVERSES
374	c=======================================================================
375	! Initialize global tracer indexes (stored in tracer.h)
376	call initracer()
377	! Load parameters from run.def file
378	CALL defrun_new( 99, .TRUE. )
379	CALL iniconst
380	CALL inigeom
381	idum=-1
382	idum=0
383
384	c Initialisation coordonnees /aires
385	c -------------------------------
386	! Note: rlatu(:) and rlonv(:) are commons defined in "comgeom.h"
387	! rlatu() and rlonv() are given in radians
388	latfi(1)=rlatu(1)
389	lonfi(1)=0.
390	DO j=2,jjm
391	DO i=1,iim
392	latfi((j-2)*iim+1+i)=rlatu(j)
393	lonfi((j-2)*iim+1+i)=rlonv(i)
394	ENDDO
395	ENDDO
396	latfi(ngridmx)=rlatu(jjp1)
397	lonfi(ngridmx)=0.
398
399	! build airefi(), mesh area on physics grid
400	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
401	! Poles are single points on physics grid
402	airefi(1)=airefi(1)*iim
403	airefi(ngridmx)=airefi(ngridmx)*iim
404
405	c=======================================================================
406	c lecture topographie, albedo, inertie thermique, relief sous-maille
407	c=======================================================================
408
409	if (choix_1.ne.1) then ! pour ne pas avoir besoin du fichier
410	! surface.dat dans le cas des start
411
412	c do while((relief(1:3).ne.'mol').AND.(relief(1:3).ne.'pla'))
413	c write(,)
414	c write(,) 'choix du relief (mola,pla)'
415	c write(,) '(Topographie MGS MOLA, plat)'
416	c read(*,fmt='(a3)') relief
417	relief="mola"
418	c enddo
419
420	CALL datareadnc(relief,phis,alb,surfith,zmeaS,zstdS,zsigS,zgamS,
421	. ztheS)
422
423	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
424	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,surfith,surfithfi)
425	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
426	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zmeaS,zmea)
427	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zstdS,zstd)
428	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zsigS,zsig)
429	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zgamS,zgam)
430	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ztheS,zthe)
431
432	endif ! of if (choix_1.ne.1)
433
434
435	c=======================================================================
436	c Lecture des fichiers (start ou start_archive)
437	c=======================================================================
438
439	if (choix_1.eq.0) then
440
441	write(,) 'Reading file START_ARCHIVE'
442	CALL lect_start_archive(date,tsurf,tsoil,emis,q2,
443	. t,ucov,vcov,ps,teta,phisold_newgrid,q,qsurf,
444	& surfith,nid)
445	write(,) "OK, read start_archive file"
446	! copy soil thermal inertia
447	ithfi(:,:)=inertiedat(:,:)
448
449	ierr= NF_CLOSE(nid)
450
451	else if (choix_1.eq.1) then
452	!do nothing, start and startfi have already been read
453	else
454	CALL exit(1)
455	endif
456
457	dtvr = daysec/FLOAT(day_step)
458	dtphys = dtvr * FLOAT(iphysiq)
459
460	c=======================================================================
461	c
462	c=======================================================================
463
464	do ! infinite loop on list of changes
465
466	write(,)
467	write(,)
468	write(,) 'List of possible changes :'
469	write(,) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
470	write(,)
471	write(,) 'flat : no topography ("aquaplanet")'
472	write(,) 'nuketharsis : no Tharsis bulge'
473	write(,) 'bilball : uniform albedo and thermal inertia'
474	write(,) 'coldspole : cold subsurface and high albedo at S.pole'
475	write(,) 'qname : change tracer name'
476	write(,) 'q=0 : ALL tracer =zero'
477	write(,) 'q=x : give a specific uniform value to one tracer'
478	write(,) 'q=profile : specify a profile for a tracer'
479	! write(,) 'ini_q : tracers initialisation for chemistry, water an
480	! $d ice '
481	! write(,) 'ini_q-H2O : tracers initialisation for chemistry and
482	! $ice '
483	! write(,) 'ini_q-iceH2O : tracers initialisation for chemistry on
484	! $ly '
485	write(,) 'noglacier : Remove tropical H2O ice if \|lat\|<45'
486	write(,) 'watercapn : H20 ice on permanent N polar cap '
487	write(,) 'watercaps : H20 ice on permanent S polar cap '
488	write(,) 'noacglac : H2O ice across Noachis Terra'
489	write(,) 'oborealis : H2O ice across Vastitas Borealis'
490	write(,) 'iceball : Thick ice layer all over surface'
491	write(,) 'wetstart : start with a wet atmosphere'
492	write(,) 'isotherm : Isothermal Temperatures, wind set to zero'
493	write(,) 'radequi : Earth-like radiative equilibrium temperature
494	$ profile (lat-alt) and winds set to zero'
495	write(,) 'coldstart : Start X K above the CO2 frost point and
496	$set wind to zero (assumes 100% CO2)'
497	write(,) 'co2ice=0 : remove CO2 polar cap'
498	write(,) 'ptot : change total pressure'
499	write(,) 'emis : change surface emissivity'
500	write(,) 'therm_ini_s : Set soil thermal inertia to reference sur
501	&face values'
502
503	write(,)
504	write(,) 'Change to perform ?'
505	write(,) ' (enter keyword or return to end)'
506	write(,)
507
508	read(*,fmt='(a20)') modif
509	if (modif(1:1) .eq. ' ') exit ! exit loop on changes
510
511	write(,)
512	write(,) trim(modif) , ' : '
513
514	c 'flat : no topography ("aquaplanet")'
515	c -------------------------------------
516	if (trim(modif) .eq. 'flat') then
517	c set topo to zero
518	z_reel(1:iip1,1:jjp1)=0
519	phis(1:iip1,1:jjp1)=g*z_reel(1:iip1,1:jjp1)
520	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
521	write(,) 'topography set to zero.'
522	write(,) 'WARNING : the subgrid topography parameters',
523	& ' were not set to zero ! => set calllott to F'
524
525	c Choice of surface pressure
526	yes=' '
527	do while ((yes.ne.'y').and.(yes.ne.'n'))
528	write(,) 'Do you wish to choose homogeneous surface',
529	& 'pressure (y) or let newstart interpolate ',
530	& ' the previous field (n)?'
531	read(*,fmt='(a)') yes
532	end do
533	if (yes.eq.'y') then
534	flagps0=.true.
535	write(,) 'New value for ps (Pa) ?'
536	201 read(,,iostat=ierr) patm
537	if(ierr.ne.0) goto 201
538	write(,)
539	write(,) ' new ps everywhere (Pa) = ', patm
540	write(,)
541	do j=1,jjp1
542	do i=1,iip1
543	ps(i,j)=patm
544	enddo
545	enddo
546	end if
547
548	c 'nuketharsis : no tharsis bulge for Early Mars'
549	c ---------------------------------------------
550	else if (trim(modif) .eq. 'nuketharsis') then
551
552	DO j=1,jjp1
553	DO i=1,iim
554	ig=1+(j-2)*iim +i
555	if(j.eq.1) ig=1
556	if(j.eq.jjp1) ig=ngridmx
557
558	fact1=(((rlonv(i)180./pi)+100)*2 +
559	& (rlatu(j)180./pi)2)/65*2
560	fact2=exp( -fact1**2.5 )
561
562	phis(i,j) = phis(i,j) - (phis(i,j)+4000.g)fact2
563
564	! if(phis(i,j).gt.2500.*g)then
565	! if(rlatu(j)*180./pi.gt.-80.)then ! avoid chopping south polar cap
566	! phis(i,j)=2500.*g
567	! endif
568	! endif
569
570	ENDDO
571	ENDDO
572	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
573
574
575	c bilball : uniform albedo, thermal inertia
576	c -----------------------------------------
577	else if (trim(modif) .eq. 'bilball') then
578	write(,) 'constante albedo and iner.therm:'
579	write(,) 'New value for albedo (ex: 0.25) ?'
580	101 read(,,iostat=ierr) alb_bb
581	if(ierr.ne.0) goto 101
582	write(,)
583	write(,) ' uniform albedo (new value):',alb_bb
584	write(,)
585
586	write(,) 'New value for thermal inertia (eg: 247) ?'
587	102 read(,,iostat=ierr) ith_bb
588	if(ierr.ne.0) goto 102
589	write(,) 'uniform thermal inertia (new value):',ith_bb
590	DO j=1,jjp1
591	DO i=1,iip1
592	alb(i,j) = alb_bb ! albedo
593	do isoil=1,nsoilmx
594	ith(i,j,isoil) = ith_bb ! thermal inertia
595	enddo
596	END DO
597	END DO
598	! CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ith,ithfi)
599	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
600	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
601
602	c coldspole : sous-sol de la calotte sud toujours froid
603	c -----------------------------------------------------
604	else if (trim(modif) .eq. 'coldspole') then
605	write(,)'new value for the subsurface temperature',
606	& ' beneath the permanent southern polar cap ? (eg: 141 K)'
607	103 read(,,iostat=ierr) tsud
608	if(ierr.ne.0) goto 103
609	write(,)
610	write(,) ' new value of the subsurface temperature:',tsud
611	c nouvelle temperature sous la calotte permanente
612	do l=2,nsoilmx
613	tsoil(ngridmx,l) = tsud
614	end do
615
616
617	write(,)'new value for the albedo',
618	& 'of the permanent southern polar cap ? (eg: 0.75)'
619	104 read(,,iostat=ierr) albsud
620	if(ierr.ne.0) goto 104
621	write(,)
622
623	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
624	c Option 1: only the albedo of the pole is modified :
625	albfi(ngridmx)=albsud
626	write(,) 'ig=',ngridmx,' albedo perennial cap ',
627	& albfi(ngridmx)
628
629	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
630	c Option 2 A haute resolution : coordonnee de la vrai calotte ~
631	c DO j=1,jjp1
632	c DO i=1,iip1
633	c ig=1+(j-2)*iim +i
634	c if(j.eq.1) ig=1
635	c if(j.eq.jjp1) ig=ngridmx
636	c if ((rlatu(j)*180./pi.lt.-84.).and.
637	c & (rlatu(j)*180./pi.gt.-91.).and.
638	c & (rlonv(i)*180./pi.gt.-91.).and.
639	c & (rlonv(i)*180./pi.lt.0.)) then
640	cc albedo de la calotte permanente fixe a albsud
641	c alb(i,j)=albsud
642	c write(,) 'lat=',rlatu(j)*180./pi,
643	c & ' lon=',rlonv(i)*180./pi
644	cc fin de la condition sur les limites de la calotte permanente
645	c end if
646	c ENDDO
647	c ENDDO
648	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
649
650	c CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
651
652
653	c ptot : Modification of the total pressure: ice + current atmosphere
654	c -------------------------------------------------------------------
655	else if (trim(modif).eq.'ptot') then
656
657	! check if we have a co2_ice surface tracer:
658	if (igcm_co2_ice.eq.0) then
659	write(,) " No surface CO2 ice !"
660	write(,) " only atmospheric pressure will be considered!"
661	endif
662	c calcul de la pression totale glace + atm actuelle
663	patm=0.
664	airetot=0.
665	pcap=0.
666	DO j=1,jjp1
667	DO i=1,iim
668	ig=1+(j-2)*iim +i
669	if(j.eq.1) ig=1
670	if(j.eq.jjp1) ig=ngridmx
671	patm = patm + ps(i,j)*aire(i,j)
672	airetot= airetot + aire(i,j)
673	if (igcm_co2_ice.ne.0) then
674	!pcap = pcap + aire(i,j)co2ice(ig)g
675	pcap = pcap + aire(i,j)qsurf(ig,igcm_co2_ice)g
676	endif
677	ENDDO
678	ENDDO
679	ptoto = pcap + patm
680
681	print*,'Current total pressure at surface (co2 ice + atm) ',
682	& ptoto/airetot
683
684	print*,'new value?'
685	read(,) ptotn
686	ptotn=ptotn*airetot
687	patmn=ptotn-pcap
688	print*,'ptoto,patm,ptotn,patmn'
689	print*,ptoto,patm,ptotn,patmn
690	print*,'Mult. factor for pressure (atm only)', patmn/patm
691	do j=1,jjp1
692	do i=1,iip1
693	ps(i,j)=ps(i,j)*patmn/patm
694	enddo
695	enddo
696
697
698
699	c Correction pour la conservation des traceurs
700	yes=' '
701	do while ((yes.ne.'y').and.(yes.ne.'n'))
702	write(,) 'Do you wish to conserve tracer total mass (y)',
703	& ' or tracer mixing ratio (n) ?'
704	read(*,fmt='(a)') yes
705	end do
706
707	if (yes.eq.'y') then
708	write(,) 'OK : conservation of tracer total mass'
709	DO iq =1, nqmx
710	DO l=1,llm
711	DO j=1,jjp1
712	DO i=1,iip1
713	q(i,j,l,iq)=q(i,j,l,iq)*patm/patmn
714	ENDDO
715	ENDDO
716	ENDDO
717	ENDDO
718	else
719	write(,) 'OK : conservation of tracer mixing ratio'
720	end if
721
722	c Correction pour la pression au niveau de la mer
723	yes=' '
724	do while ((yes.ne.'y').and.(yes.ne.'n'))
725	write(,) 'Do you wish fix pressure at sea level (y)',
726	& ' or not (n) ?'
727	read(*,fmt='(a)') yes
728	end do
729
730	if (yes.eq.'y') then
731	write(,) 'Value?'
732	read(,,iostat=ierr) psea
733	DO i=1,iip1
734	DO j=1,jjp1
735	ps(i,j)=psea
736
737	ENDDO
738	ENDDO
739	write(,) 'psea=',psea
740	else
741	write(,) 'no'
742	end if
743
744
745	c emis : change surface emissivity (added by RW)
746	c ----------------------------------------------
747	else if (trim(modif).eq.'emis') then
748
749	print*,'new value?'
750	read(,) emisread
751
752	do i=1,ngridmx
753	emis(i)=emisread
754	enddo
755
756	c qname : change tracer name
757	c --------------------------
758	else if (trim(modif).eq.'qname') then
759	yes='y'
760	do while (yes.eq.'y')
761	write(,) 'Which tracer name do you want to change ?'
762	do iq=1,nqmx
763	write(*,'(i3,a3,a20)')iq,' : ',trim(tnom(iq))
764	enddo
765	write(*,'(a35,i3)')
766	& '(enter tracer number; between 1 and ',nqmx
767	write(,)' or any other value to quit this option)'
768	read(,) iq
769	if ((iq.ge.1).and.(iq.le.nqmx)) then
770	write(,)'Change tracer name ',trim(tnom(iq)),' to ?'
771	read(,) txt
772	tnom(iq)=txt
773	write(,)'Do you want to change another tracer name (y/n)?'
774	read(*,'(a)') yes
775	else
776	! inapropiate value of iq; quit this option
777	yes='n'
778	endif ! of if ((iq.ge.1).and.(iq.le.nqmx))
779	enddo ! of do while (yes.ne.'y')
780
781	c q=0 : set tracers to zero
782	c -------------------------
783	else if (trim(modif).eq.'q=0') then
784	c mise a 0 des q (traceurs)
785	write(,) 'Tracers set to 0 (1.E-30 in fact)'
786	DO iq =1, nqmx
787	DO l=1,llm
788	DO j=1,jjp1
789	DO i=1,iip1
790	q(i,j,l,iq)=1.e-30
791	ENDDO
792	ENDDO
793	ENDDO
794	ENDDO
795
796	c set surface tracers to zero
797	DO iq =1, nqmx
798	DO ig=1,ngridmx
799	qsurf(ig,iq)=0.
800	ENDDO
801	ENDDO
802
803	c q=x : initialise tracer manually
804	c --------------------------------
805	else if (trim(modif).eq.'q=x') then
806	write(,) 'Which tracer do you want to modify ?'
807	do iq=1,nqmx
808	write(,)iq,' : ',trim(tnom(iq))
809	enddo
810	write(,) '(choose between 1 and ',nqmx,')'
811	read(,) iq
812	write(,)'mixing ratio of tracer ',trim(tnom(iq)),
813	& ' ? (kg/kg)'
814	read(,) val
815	DO l=1,llm
816	DO j=1,jjp1
817	DO i=1,iip1
818	q(i,j,l,iq)=val
819	ENDDO
820	ENDDO
821	ENDDO
822	write(,) 'SURFACE value of tracer ',trim(tnom(iq)),
823	& ' ? (kg/m2)'
824	read(,) val
825	DO ig=1,ngridmx
826	qsurf(ig,iq)=val
827	ENDDO
828
829	c q=profile : initialize tracer with a given profile
830	c --------------------------------------------------
831	else if (trim(modif) .eq. 'q=profile') then
832	write(,) 'Tracer profile will be sought in ASCII file'
833	write(,) "'profile_tracer' where 'tracer' is tracer name"
834	write(,) "(one value per line in file; starting with"
835	write(,) "surface value, the 1st atmospheric layer"
836	write(,) "followed by 2nd, etc. up to top of atmosphere)"
837	write(,) 'Which tracer do you want to set?'
838	do iq=1,nqmx
839	write(,)iq,' : ',trim(tnom(iq))
840	enddo
841	write(,) '(choose between 1 and ',nqmx,')'
842	read(,) iq
843	if ((iq.lt.1).or.(iq.gt.nqmx)) then
844	! wrong value for iq, go back to menu
845	write(,) "wrong input value:",iq
846	cycle
847	endif
848	! look for input file 'profile_tracer'
849	txt="profile_"//trim(tnom(iq))
850	open(41,file=trim(txt),status='old',form='formatted',
851	& iostat=ierr)
852	if (ierr.eq.0) then
853	! OK, found file 'profile_...', load the profile
854	do l=1,llm+1
855	read(41,*,iostat=ierr) profile(l)
856	if (ierr.ne.0) then ! something went wrong
857	exit ! quit loop
858	endif
859	enddo
860	if (ierr.eq.0) then
861	! initialize tracer values
862	qsurf(:,iq)=profile(1)
863	do l=1,llm
864	q(:,:,l,iq)=profile(l+1)
865	enddo
866	write(,)'OK, tracer ',trim(tnom(iq)),' initialized ',
867	& 'using values from file ',trim(txt)
868	else
869	write(,)'problem reading file ',trim(txt),' !'
870	write(,)'No modifications to tracer ',trim(tnom(iq))
871	endif
872	else
873	write(,)'Could not find file ',trim(txt),' !'
874	write(,)'No modifications to tracer ',trim(tnom(iq))
875	endif
876
877
878	c wetstart : wet atmosphere with a north to south gradient
879	c --------------------------------------------------------
880	else if (trim(modif) .eq. 'wetstart') then
881	! check that there is indeed a water vapor tracer
882	if (igcm_h2o_vap.eq.0) then
883	write(,) "No water vapour tracer! Can't use this option"
884	stop
885	endif
886	DO l=1,llm
887	DO j=1,jjp1
888	DO i=1,iip1
889	q(i,j,l,igcm_h2o_vap)=150.e-6 * (rlatu(j)+pi/2.) / pi
890	ENDDO
891	ENDDO
892	ENDDO
893
894	write(,) 'Water mass mixing ratio at north pole='
895	* ,q(1,1,1,igcm_h2o_vap)
896	write(,) '---------------------------south pole='
897	* ,q(1,jjp1,1,igcm_h2o_vap)
898
899	c noglacier : remove tropical water ice (to initialize high res sim)
900	c --------------------------------------------------
901	else if (trim(modif) .eq. 'noglacier') then
902	if (igcm_h2o_ice.eq.0) then
903	write(,) "No water ice tracer! Can't use this option"
904	stop
905	endif
906	do ig=1,ngridmx
907	j=(ig-2)/iim +2
908	if(ig.eq.1) j=1
909	write(,) 'OK: remove surface ice for \|lat\|<45'
910	if (abs(rlatu(j)*180./pi).lt.45.) then
911	qsurf(ig,igcm_h2o_ice)=0.
912	end if
913	end do
914
915
916	c watercapn : H20 ice on permanent northern cap
917	c --------------------------------------------------
918	else if (trim(modif) .eq. 'watercapn') then
919	if (igcm_h2o_ice.eq.0) then
920	write(,) "No water ice tracer! Can't use this option"
921	stop
922	endif
923
924	DO j=1,jjp1
925	DO i=1,iim
926	ig=1+(j-2)*iim +i
927	if(j.eq.1) ig=1
928	if(j.eq.jjp1) ig=ngridmx
929
930	if (rlatu(j)*180./pi.gt.80.) then
931	qsurf(ig,igcm_h2o_ice)=3.4e3
932	!do isoil=1,nsoilmx
933	! ith(i,j,isoil) = 18000. ! thermal inertia
934	!enddo
935	write(,)' ==> Ice mesh North boundary (deg)= ',
936	& rlatv(j-1)*180./pi
937	end if
938	ENDDO
939	ENDDO
940	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
941
942	c$$$ do ig=1,ngridmx
943	c$$$ j=(ig-2)/iim +2
944	c$$$ if(ig.eq.1) j=1
945	c$$$ if (rlatu(j)*180./pi.gt.80.) then
946	c$$$
947	c$$$ qsurf(ig,igcm_h2o_ice)=1.e5
948	c$$$ qsurf(ig,igcm_h2o_vap)=0.0!1.e5
949	c$$$
950	c$$$ write(,) 'ig=',ig,' H2O ice mass (kg/m2)= ',
951	c$$$ & qsurf(ig,igcm_h2o_ice)
952	c$$$
953	c$$$ write(,)' ==> Ice mesh South boundary (deg)= ',
954	c$$$ & rlatv(j)*180./pi
955	c$$$ end if
956	c$$$ enddo
957
958	c watercaps : H20 ice on permanent southern cap
959	c -------------------------------------------------
960	else if (trim(modif) .eq. 'watercaps') then
961	if (igcm_h2o_ice.eq.0) then
962	write(,) "No water ice tracer! Can't use this option"
963	stop
964	endif
965
966	DO j=1,jjp1
967	DO i=1,iim
968	ig=1+(j-2)*iim +i
969	if(j.eq.1) ig=1
970	if(j.eq.jjp1) ig=ngridmx
971
972	if (rlatu(j)*180./pi.lt.-80.) then
973	qsurf(ig,igcm_h2o_ice)=3.4e3
974	!do isoil=1,nsoilmx
975	! ith(i,j,isoil) = 18000. ! thermal inertia
976	!enddo
977	write(,)' ==> Ice mesh South boundary (deg)= ',
978	& rlatv(j-1)*180./pi
979	end if
980	ENDDO
981	ENDDO
982	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
983
984	c$$$ do ig=1,ngridmx
985	c$$$ j=(ig-2)/iim +2
986	c$$$ if(ig.eq.1) j=1
987	c$$$ if (rlatu(j)*180./pi.lt.-80.) then
988	c$$$ qsurf(ig,igcm_h2o_ice)=1.e5
989	c$$$ qsurf(ig,igcm_h2o_vap)=0.0 !1.e5
990	c$$$
991	c$$$ write(,) 'ig=',ig,' H2O ice mass (kg/m2)= ',
992	c$$$ & qsurf(ig,igcm_h2o_ice)
993	c$$$ write(,)' ==> Ice mesh North boundary (deg)= ',
994	c$$$ & rlatv(j-1)*180./pi
995	c$$$ end if
996	c$$$ enddo
997
998
999	c noacglac : H2O ice across highest terrain
1000	c --------------------------------------------
1001	else if (trim(modif) .eq. 'noacglac') then
1002	if (igcm_h2o_ice.eq.0) then
1003	write(,) "No water ice tracer! Can't use this option"
1004	stop
1005	endif
1006	DO j=1,jjp1
1007	DO i=1,iim
1008	ig=1+(j-2)*iim +i
1009	if(j.eq.1) ig=1
1010	if(j.eq.jjp1) ig=ngridmx
1011
1012	if(phis(i,j).gt.1000.*g)then
1013	alb(i,j) = 0.5 ! snow value
1014	do isoil=1,nsoilmx
1015	ith(i,j,isoil) = 18000. ! thermal inertia
1016	! this leads to rnat set to 'ocean' in physiq.F90
1017	! actually no, because it is soil not surface
1018	enddo
1019	endif
1020	ENDDO
1021	ENDDO
1022	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
1023	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
1024	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
1025
1026
1027
1028	c oborealis : H2O oceans across Vastitas Borealis
1029	c -----------------------------------------------
1030	else if (trim(modif) .eq. 'oborealis') then
1031	if (igcm_h2o_ice.eq.0) then
1032	write(,) "No water ice tracer! Can't use this option"
1033	stop
1034	endif
1035	DO j=1,jjp1
1036	DO i=1,iim
1037	ig=1+(j-2)*iim +i
1038	if(j.eq.1) ig=1
1039	if(j.eq.jjp1) ig=ngridmx
1040
1041	if(phis(i,j).lt.-4000.*g)then
1042	! if( (phis(i,j).lt.-4000.*g)
1043	! & .and. (rlatu(j)*180./pi.lt.0.) )then ! south hemisphere only
1044	! phis(i,j)=-8000.0*g ! proper ocean
1045	phis(i,j)=-4000.0*g ! scanty ocean
1046
1047	alb(i,j) = 0.07 ! oceanic value
1048	do isoil=1,nsoilmx
1049	ith(i,j,isoil) = 18000. ! thermal inertia
1050	! this leads to rnat set to 'ocean' in physiq.F90
1051	! actually no, because it is soil not surface
1052	enddo
1053	endif
1054	ENDDO
1055	ENDDO
1056	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
1057	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
1058	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
1059
1060	c iborealis : H2O ice in Northern plains
1061	c --------------------------------------
1062	else if (trim(modif) .eq. 'iborealis') then
1063	if (igcm_h2o_ice.eq.0) then
1064	write(,) "No water ice tracer! Can't use this option"
1065	stop
1066	endif
1067	DO j=1,jjp1
1068	DO i=1,iim
1069	ig=1+(j-2)*iim +i
1070	if(j.eq.1) ig=1
1071	if(j.eq.jjp1) ig=ngridmx
1072
1073	if(phis(i,j).lt.-4000.*g)then
1074	!qsurf(ig,igcm_h2o_ice)=1.e3
1075	qsurf(ig,igcm_h2o_ice)=241.4 ! to make total 33 kg m^-2
1076	endif
1077	ENDDO
1078	ENDDO
1079	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
1080	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
1081	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
1082
1083
1084	c oceanball : H2O liquid everywhere
1085	c ----------------------------
1086	else if (trim(modif) .eq. 'oceanball') 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	qsurf(ig,igcm_h2o_vap)=0.0 ! for ocean, this is infinite source
1098	qsurf(ig,igcm_h2o_ice)=0.0
1099	alb(i,j) = 0.07 ! ocean value
1100
1101	do isoil=1,nsoilmx
1102	ith(i,j,isoil) = 18000. ! thermal inertia
1103	!ith(i,j,isoil) = 50. ! extremely low for test
1104	! this leads to rnat set to 'ocean' in physiq.F90
1105	enddo
1106
1107	ENDDO
1108	ENDDO
1109	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
1110	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
1111	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
1112
1113	c iceball : H2O ice everywhere
1114	c ----------------------------
1115	else if (trim(modif) .eq. 'iceball') then
1116	if (igcm_h2o_ice.eq.0) then
1117	write(,) "No water ice tracer! Can't use this option"
1118	stop
1119	endif
1120	DO j=1,jjp1
1121	DO i=1,iim
1122	ig=1+(j-2)*iim +i
1123	if(j.eq.1) ig=1
1124	if(j.eq.jjp1) ig=ngridmx
1125
1126	qsurf(ig,igcm_h2o_vap)=-50. ! for ocean, this is infinite source
1127	qsurf(ig,igcm_h2o_ice)=50. ! == to 0.05 m of oceanic ice
1128	alb(i,j) = 0.6 ! ice albedo value
1129
1130	do isoil=1,nsoilmx
1131	!ith(i,j,isoil) = 18000. ! thermal inertia
1132	! this leads to rnat set to 'ocean' in physiq.F90
1133	enddo
1134
1135	ENDDO
1136	ENDDO
1137	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
1138	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
1139
1140	c isotherm : Isothermal temperatures and no winds
1141	c -----------------------------------------------
1142	else if (trim(modif) .eq. 'isotherm') then
1143
1144	write(,)'Isothermal temperature of the atmosphere,
1145	& surface and subsurface'
1146	write(,) 'Value of this temperature ? :'
1147	203 read(,,iostat=ierr) Tiso
1148	if(ierr.ne.0) goto 203
1149
1150	tsurf(1:ngridmx)=Tiso
1151
1152	tsoil(1:ngridmx,1:nsoilmx)=Tiso
1153
1154	Tset(1:iip1,1:jjp1,1:llm)=Tiso
1155	flagtset=.true.
1156
1157	ucov(1:iip1,1:jjp1,1:llm)=0
1158	vcov(1:iip1,1:jjm,1:llm)=0
1159	q2(1:ngridmx,1:nlayermx+1)=0
1160
1161	c radequi : Radiative equilibrium profile of temperatures and no winds
1162	c --------------------------------------------------------------------
1163	else if (trim(modif) .eq. 'radequi') then
1164
1165	write(,)'radiative equilibrium temperature profile'
1166
1167	do ig=1, ngridmx
1168	teque= 335.0-60.0sin(latfi(ig))sin(latfi(ig))-
1169	& 10.0cos(latfi(ig))cos(latfi(ig))
1170	tsurf(ig) = MAX(220.0,teque)
1171	end do
1172	do l=2,nsoilmx
1173	do ig=1, ngridmx
1174	tsoil(ig,l) = tsurf(ig)
1175	end do
1176	end do
1177	DO j=1,jjp1
1178	DO i=1,iim
1179	Do l=1,llm
1180	teque=335.-60.0sin(rlatu(j))sin(rlatu(j))-
1181	& 10.0cos(rlatu(j))cos(rlatu(j))
1182	Tset(i,j,l)=MAX(220.0,teque)
1183	end do
1184	end do
1185	end do
1186	flagtset=.true.
1187	ucov(1:iip1,1:jjp1,1:llm)=0
1188	vcov(1:iip1,1:jjm,1:llm)=0
1189	q2(1:ngridmx,1:nlayermx+1)=0
1190
1191	c coldstart : T set 1K above CO2 frost point and no winds
1192	c ------------------------------------------------
1193	else if (trim(modif) .eq. 'coldstart') then
1194
1195	write(,)'set temperature of the atmosphere,'
1196	&,'surface and subsurface how many degrees above CO2 frost point?'
1197	204 read(,,iostat=ierr) Tabove
1198	if(ierr.ne.0) goto 204
1199
1200	DO j=1,jjp1
1201	DO i=1,iim
1202	ig=1+(j-2)*iim +i
1203	if(j.eq.1) ig=1
1204	if(j.eq.jjp1) ig=ngridmx
1205	tsurf(ig) = (-3167.8)/(log(.01*ps(i,j))-23.23)+Tabove
1206	END DO
1207	END DO
1208	do l=1,nsoilmx
1209	do ig=1, ngridmx
1210	tsoil(ig,l) = tsurf(ig)
1211	end do
1212	end do
1213	DO j=1,jjp1
1214	DO i=1,iim
1215	Do l=1,llm
1216	pp = aps(l) +bps(l)*ps(i,j)
1217	Tset(i,j,l)=(-3167.8)/(log(.01*pp)-23.23)+Tabove
1218	end do
1219	end do
1220	end do
1221
1222	flagtset=.true.
1223	ucov(1:iip1,1:jjp1,1:llm)=0
1224	vcov(1:iip1,1:jjm,1:llm)=0
1225	q2(1:ngridmx,1:nlayermx+1)=0
1226
1227
1228	c co2ice=0 : remove CO2 polar ice caps'
1229	c ------------------------------------------------
1230	else if (trim(modif) .eq. 'co2ice=0') then
1231	! check that there is indeed a co2_ice tracer ...
1232	if (igcm_co2_ice.ne.0) then
1233	do ig=1,ngridmx
1234	!co2ice(ig)=0
1235	qsurf(ig,igcm_co2_ice)=0
1236	emis(ig)=emis(ngridmx/2)
1237	end do
1238	else
1239	write(,) "Can't remove CO2 ice!! (no co2_ice tracer)"
1240	endif
1241
1242	! therm_ini_s: (re)-set soil thermal inertia to reference surface values
1243	! ----------------------------------------------------------------------
1244
1245	else if (trim(modif).eq.'therm_ini_s') then
1246	! write(,)"surfithfi(1):",surfithfi(1)
1247	do isoil=1,nsoilmx
1248	inertiedat(1:ngridmx,isoil)=surfithfi(1:ngridmx)
1249	enddo
1250	write(,)'OK: Soil thermal inertia has been reset to referenc
1251	&e surface values'
1252	! write(,)"inertiedat(1,1):",inertiedat(1,1)
1253	ithfi(:,:)=inertiedat(:,:)
1254	! recast ithfi() onto ith()
1255	call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
1256	! Check:
1257	! do i=1,iip1
1258	! do j=1,jjp1
1259	! do isoil=1,nsoilmx
1260	! write(77,*) i,j,isoil," ",ith(i,j,isoil)
1261	! enddo
1262	! enddo
1263	! enddo
1264
1265
1266	else
1267	write(,) ' Unknown (misspelled?) option!!!'
1268	end if ! of if (trim(modif) .eq. '...') elseif ...
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280	enddo ! of do ! infinite loop on liste of changes
1281
1282	999 continue
1283
1284
1285	c=======================================================================
1286	c Initialisation for cloud fraction and oceanic ice (added by BC 2010)
1287	c=======================================================================
1288	DO ig=1, ngridmx
1289	totalfrac(ig)=0.5
1290	DO l=1,llm
1291	cloudfrac(ig,l)=0.5
1292	ENDDO
1293	! Ehouarn, march 2012: also add some initialisation for hice
1294	hice(ig)=0.0
1295	ENDDO
1296
1297	c========================================================
1298
1299	! DO ig=1,ngridmx
1300	! IF(tsurf(ig) .LT. 273.16-1.8) THEN
1301	! hice(ig)=(273.16-1.8-tsurf(ig))/(273.16-1.8-240)*1.
1302	! hice(ig)=min(hice(ig),1.0)
1303	! ENDIF
1304	! ENDDO
1305
1306	c=======================================================================
1307	c Correct pressure on the new grid (menu 0)
1308	c=======================================================================
1309
1310
1311	if ((choix_1.eq.0).and.(.not.flagps0)) then
1312	r = 1000.*8.31/mugaz
1313
1314	do j=1,jjp1
1315	do i=1,iip1
1316	ps(i,j) = ps(i,j) *
1317	. exp((phisold_newgrid(i,j)-phis(i,j)) /
1318	. (t(i,j,1) * r))
1319	end do
1320	end do
1321
1322	c periodicite de ps en longitude
1323	do j=1,jjp1
1324	ps(1,j) = ps(iip1,j)
1325	end do
1326	end if
1327
1328
1329	c=======================================================================
1330	c=======================================================================
1331
1332	c=======================================================================
1333	c Initialisation de la physique / ecriture de newstartfi :
1334	c=======================================================================
1335
1336
1337	CALL inifilr
1338	CALL pression(ip1jmp1, ap, bp, ps, p3d)
1339
1340	c-----------------------------------------------------------------------
1341	c Initialisation pks:
1342	c-----------------------------------------------------------------------
1343
1344	CALL exner_hyb(ip1jmp1, ps, p3d, beta, pks, pk, pkf)
1345	! Calcul de la temperature potentielle teta
1346
1347	if (flagtset) then
1348	DO l=1,llm
1349	DO j=1,jjp1
1350	DO i=1,iim
1351	teta(i,j,l) = Tset(i,j,l) * cpp/pk(i,j,l)
1352	ENDDO
1353	teta (iip1,j,l)= teta (1,j,l)
1354	ENDDO
1355	ENDDO
1356	else if (choix_1.eq.0) then
1357	DO l=1,llm
1358	DO j=1,jjp1
1359	DO i=1,iim
1360	teta(i,j,l) = t(i,j,l) * cpp/pk(i,j,l)
1361	ENDDO
1362	teta (iip1,j,l)= teta (1,j,l)
1363	ENDDO
1364	ENDDO
1365	endif
1366
1367	C Calcul intermediaire
1368	c
1369	if (choix_1.eq.0) then
1370	CALL massdair( p3d, masse )
1371	c
1372	print *,' ALPHAX ',alphax
1373
1374	DO l = 1, llm
1375	DO i = 1, iim
1376	xppn(i) = aire( i, 1 ) * masse( i , 1 , l )
1377	xpps(i) = aire( i,jjp1 ) * masse( i , jjp1 , l )
1378	ENDDO
1379	xpn = SUM(xppn)/apoln
1380	xps = SUM(xpps)/apols
1381	DO i = 1, iip1
1382	masse( i , 1 , l ) = xpn
1383	masse( i , jjp1 , l ) = xps
1384	ENDDO
1385	ENDDO
1386	endif
1387	phis(iip1,:) = phis(1,:)
1388
1389	CALL inidissip ( lstardis, nitergdiv, nitergrot, niterh,
1390	* tetagdiv, tetagrot , tetatemp )
1391	itau=0
1392	if (choix_1.eq.0) then
1393	day_ini=int(date)
1394	endif
1395	c
1396	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
1397
1398	CALL caldyn0( itau,ucov,vcov,teta,ps,masse,pk,phis ,
1399	* phi,w, pbaru,pbarv,day_ini+time )
1400
1401
1402	CALL dynredem0("restart.nc",day_ini,anneeref,phis,nqmx)
1403	CALL dynredem1("restart.nc",0.0,vcov,ucov,teta,q,nqmx,masse,ps)
1404	C
1405	C Ecriture etat initial physique
1406	C
1407
1408	! do ig=1,ngridmx
1409	! print*,'surface ice in newstart=',qsurf(ig,igcm_h2o_ice)
1410	! print*,'surface liq in newstart=',qsurf(ig,igcm_h2o_vap)
1411	! enddo
1412	! stop
1413
1414
1415	call physdem1("restartfi.nc",lonfi,latfi,nsoilmx,nqmx,
1416	. dtphys,float(day_ini),
1417	. time,tsurf,tsoil,emis,q2,qsurf,
1418	. airefi,albfi,ithfi,zmea,zstd,zsig,zgam,zthe,
1419	. cloudfrac,totalfrac,hice)
1420
1421	c=======================================================================
1422	c Formats
1423	c=======================================================================
1424
1425	1 FORMAT(//10x,'la valeur de im =',i4,2x,'lue sur le fichier de dema
1426	*rrage est differente de la valeur parametree iim =',i4//)
1427	2 FORMAT(//10x,'la valeur de jm =',i4,2x,'lue sur le fichier de dema
1428	*rrage est differente de la valeur parametree jjm =',i4//)
1429	3 FORMAT(//10x,'la valeur de lllm =',i4,2x,'lue sur le fichier demar
1430	*rage est differente de la valeur parametree llm =',i4//)
1431
1432	end
1433

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

Download in other formats:

Original Format