Context Navigation

newstart.F @ 1711

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

Mars GCM

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

File size: 62.1 KB

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

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

Download in other formats:

Original Format