Context Navigation

newstart.F @ 1241

Last change on this file since 1241 was 1241, checked in by aslmd, 11 years ago
LMDZ.MARS. bug introduced in rev1233 in newstart: ngrid is not present in newstart. now fixed.
File size: 62.3 KB

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

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

Download in other formats:

Original Format