Context Navigation

newstart.F @ 224

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

Mars GCM:

Implemented using 'z0' roughness length map (important: 'z0' reference

field is in datafile surface.nc, which has also been updated).

made z0 a z0(ngridmx) array and moved 'z0' from 'planete.h' to 'surfdat.h'; added a 'z0_default' (common in surfdat.h) corresponding to the 'control' array value (contole(19) in startfi.nc).
adapted 'tabfi.F' to use 'z0_default'.
adapted 'phyetat0.F' to look for a 'z0' field in startfi.nc. If 'z0' is not found in the startfi.nc file, then the uniform default value (z0_default) is used.
modified 'physdem1.F' to write 'z0' field to restart.nc
adapted use of z0() in 'physiq.F' (diagnostic computation of surface stress), 'vdifc.F' and 'vdif_cd.F'.
adapted 'dustdevil.F' to use 'z0_default'.
'testphys1d.F' now uses 'z0_default', and the value to use can be set in run.def (with "z0=TheValueYouWant?").
modified 'datareadnc.F' to load reference map of 'z0' from surface.nc, and added a 'z0' option in 'newstart.F' to force a uniform value of z0. Note that the use of the z0 map is automatic when using newstart, but only when it loads a start_archive.nc file.

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

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

Download in other formats:

Original Format