Context Navigation

newstart.old2 @ 205

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

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

Download in other formats:

Original Format