Context Navigation

newstart.F @ 165

Last change on this file since 165 was 164, checked in by emillour, 14 years ago

Mars GCM:

Updates and corrections (to enable compiling/running in debug mode with

ifort)

removed option "-free-form" from makegcm_ifort and set mod_loc_dir="." so that module files (produced in local directory by ifort) are moved to LIBO
updated initracer.F, physdem1.F, physiq.F, inichim_newstart.F to avoid referencing out-of-bound array indexes (even if unused)
cosmetic updates on inwrite.F, datareadnc.F
updated newstart.F to initialize and use 'datadir' when looking for files
corrected bug on interpolation of sub-surface temperatures in lect_start_archive.F

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

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

Download in other formats:

Original Format