Context Navigation

newstart.F @ 1421

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

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

Download in other formats:

Original Format