Context Navigation

newstart.F @ 1422

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

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

Download in other formats:

Original Format