Context Navigation

newstart.F @ 1711

Last change on this file since 1711 was 1711, checked in by mvals, 8 years ago

Mars GCM

implementation of a sub-grid water cloud fraction scheme (by A. Pottier): If CLFvarying is set to true in callphys.def, then the sub-grid cloud fraction routine is applied in watercloud.F and aeropacity.F.
accordingly modified files are: watercloud.F, aeropacity.F, callradite.F, conf_phys.F, phyetat0.F, phyredem.F90, physiq_mod.F, testphys1d.F, callkeys.h, newtsart.F, start2archive.F, lect_start_archive.F
added file: tcondwater.F90, used by watercloud.F to calculate the condensation temperature of water
watercloud.F, aeropacity.F, callradite.F are converted to module files as watercloud_mod.F, aeropacity_mod.F, callradite_mod.F

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

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

Download in other formats:

Original Format