Context Navigation

newstart.F @ 1231

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

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

Download in other formats:

Original Format