Context Navigation

newstart.F @ 586

Last change on this file since 586 was 535, checked in by rwordsworth, 13 years ago
fixed kappa bug
File size: 65.8 KB

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

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

Download in other formats:

Original Format