Context Navigation

newstart.F @ 1564

Last change on this file since 1564 was 1543, checked in by emillour, 9 years ago

All models: Further adaptations to keep up with changes in LMDZ5 concerning
physics/dynamics separation:

dyn3d:
adapted gcm.F so that all physics initializations are now done in iniphysiq.

dyn3dpar:
adapted gcm.F so that all physics initializations are now done in iniphysiq.
updated calfis_p.F to follow up with changes.
copied over updated "bands.F90" from LMDZ5.

dynphy_lonlat:
calfis_p.F90, mod_interface_dyn_phys.F90, follow up of changes in phy_common/mod_* routines

phy_common:
added "geometry_mod.F90" to store information about the grid (replaces phy*/comgeomphy.F90) and give variables friendlier names: rlond => longitude , rlatd => latitude, airephy => cell_area, cuphy => dx , cvphy => dy
added "physics_distribution_mod.F90"
updated "mod_grid_phy_lmdz.F90", "mod_phys_lmdz_mpi_data.F90", "mod_phys_lmdz_para.F90", "mod_phys_lmdz_mpi_transfert.F90", "mod_grid_phy_lmdz.F90", "mod_phys_lmdz_omp_data.F90", "mod_phys_lmdz_omp_transfert.F90", "write_field_phy.F90" and "ioipsl_getin_p_mod.F90" to LMDZ5 versions.

phy[venus/titan/mars/std]:
removed "init_phys_lmdz.F90", "comgeomphy.F90"; adapted routines to use geometry_mod (longitude, latitude, cell_area, etc.)

File size: 61.5 KB

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

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

Download in other formats:

Original Format