Context Navigation

newstart.F @ 1375

Last change on this file since 1375 was 1266, checked in by aslmd, 11 years ago

LMDZ.MARS
IMPORTANT CHANGE

Remove all reference/use of nlayermx and dimphys.h
Made use of automatic arrays whenever arrays are needed with dimension nlayer
Remove lots of obsolete reference to dimensions.h
Converted iono.h and param_v4.h into corresponding modules

(with embedded subroutine to allocate arrays)
(no arrays allocated if thermosphere not used)

Deleted param.h and put contents into module param_v4_h
Adapted testphys1d, newstart, etc...
Made DATA arrays in param_read to be initialized by subroutine

fill_data_thermos in module param_v4_h

Optimized computations in paramfoto_compact (twice less dlog10 calculations)
Checked consistency before/after modification in debug mode
Checked performance is not impacted (same as before)

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

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

Download in other formats:

Original Format