Context Navigation

newstart.F @ 1789

Last change on this file since 1789 was 1789, checked in by jvatant, 7 years ago
Added the surface methane tank and put it in start files --JVO
File size: 47.0 KB

Line
1	C======================================================================
2	PROGRAM newstart
3	c=======================================================================
4	c
5	c
6	c Auteur: Christophe Hourdin/Francois Forget/Yann Wanherdrick
7	c ------
8	c Derniere modif : 12/03
9	c
10	c
11	c Objet: Create or modify the initial state for the LMD Mars GCM
12	c ----- (fichiers NetCDF start et startfi)
13	c
14	c
15	c=======================================================================
16
17	use mod_phys_lmdz_para, only: is_parallel, is_sequential,
18	& is_mpi_root, is_omp_root,
19	& is_master
20	use infotrac, only: infotrac_init, nqtot, tname
21	USE comsoil_h, ONLY: nsoilmx, layer, mlayer, inertiedat
22	USE surfdat_h, ONLY: phisfi, albedodat,
23	& zmea, zstd, zsig, zgam, zthe
24	use datafile_mod, only: datadir, surfdir
25	use ioipsl_getin_p_mod, only: getin_p
26	use control_mod, only: day_step, iphysiq, anneeref, planet_type
27	use phyredem, only: physdem0, physdem1
28	use iostart, only: open_startphy
29	use filtreg_mod, only: inifilr
30	USE mod_const_mpi, ONLY: COMM_LMDZ
31	USE comvert_mod, ONLY: ap,bp,aps,bps,pa,preff
32	USE comconst_mod, ONLY: lllm,daysec,dtvr,dtphys,cpp,kappa,
33	. rad,omeg,g,r,pi
34	USE serre_mod, ONLY: alphax
35	USE temps_mod, ONLY: day_ini
36	USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
37	use tabfi_mod, only: tabfi
38	use iniphysiq_mod, only: iniphysiq
39	use phyetat0_mod, only: phyetat0
40	implicit none
41
42	include "dimensions.h"
43	integer, parameter :: ngridmx = (2+(jjm-1)*iim - 1/jjm)
44	include "paramet.h"
45	include "comgeom2.h"
46	include "comdissnew.h"
47	include "netcdf.inc"
48
49	c=======================================================================
50	c Declarations
51	c=======================================================================
52
53	c Variables dimension du fichier "start_archive"
54	c------------------------------------
55	CHARACTER relief*3
56
57
58	c Variables pour les lectures NetCDF des fichiers "start_archive"
59	c--------------------------------------------------
60	INTEGER nid_dyn, nid_fi,nid,nvarid
61	INTEGER length
62	parameter (length = 100)
63	INTEGER tab0
64	INTEGER NB_ETATMAX
65	parameter (NB_ETATMAX = 100)
66
67	REAL date
68	REAL p_rad,p_omeg,p_g,p_cpp,p_mugaz,p_daysec
69
70	c Variable histoire
71	c------------------
72	REAL vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm) ! vents covariants
73	REAL phis(iip1,jjp1)
74	REAL,ALLOCATABLE :: q(:,:,:,:) ! champs advectes
75
76	c autre variables dynamique nouvelle grille
77	c------------------------------------------
78	REAL pks(iip1,jjp1)
79	REAL w(iip1,jjp1,llm+1)
80	REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
81	! REAL dv(ip1jm,llm),du(ip1jmp1,llm)
82	! REAL dh(ip1jmp1,llm),dp(ip1jmp1)
83	REAL phi(iip1,jjp1,llm)
84
85	integer klatdat,klongdat
86	PARAMETER (klatdat=180,klongdat=360)
87
88	c Physique sur grille scalaire
89	c----------------------------
90	real zmeaS(iip1,jjp1),zstdS(iip1,jjp1)
91	real zsigS(iip1,jjp1),zgamS(iip1,jjp1),ztheS(iip1,jjp1)
92
93	c variable physique
94	c------------------
95	REAL tsurf(ngridmx) ! surface temperature
96	REAL,ALLOCATABLE :: tsoil(:,:) ! soil temperature
97	REAL emis(ngridmx) ! surface emissivity
98	real emisread ! added by RW
99	REAL,ALLOCATABLE :: qsurf(:,:)
100	REAL q2(ngridmx,llm+1)
101	! REAL rnaturfi(ngridmx)
102	real alb(iip1,jjp1),albfi(ngridmx) ! albedos
103	real,ALLOCATABLE :: ith(:,:,:),ithfi(:,:) ! thermal inertia (3D)
104	real surfith(iip1,jjp1),surfithfi(ngridmx) ! surface thermal inertia (2D)
105	REAL latfi(ngridmx),lonfi(ngridmx),airefi(ngridmx)
106
107	INTEGER i,j,l,isoil,ig,idum
108	real mugaz ! molar mass of the atmosphere
109
110	integer ierr
111
112	c Variables on the new grid along scalar points
113	c------------------------------------------------------
114	! REAL p(iip1,jjp1)
115	REAL t(iip1,jjp1,llm)
116	REAL tset(iip1,jjp1,llm)
117	real phisold_newgrid(iip1,jjp1)
118	REAL :: teta(iip1, jjp1, llm)
119	REAL :: pk(iip1,jjp1,llm)
120	REAL :: pkf(iip1,jjp1,llm)
121	REAL :: ps(iip1, jjp1)
122	REAL :: masse(iip1,jjp1,llm)
123	REAL :: xpn,xps,xppn(iim),xpps(iim)
124	REAL :: p3d(iip1, jjp1, llm+1)
125	REAL :: beta(iip1,jjp1,llm)
126	! REAL dteta(ip1jmp1,llm)
127
128	c Variable de l'ancienne grille
129	c------------------------------
130	real time
131	real tab_cntrl(100)
132	real tab_cntrl_bis(100)
133
134	c variables diverses
135	c-------------------
136	real choix_1,pp
137	character*80 fichnom
138	character*250 filestring
139	integer Lmodif,iq
140	character modif*20
141	real z_reel(iip1,jjp1)
142	real tsud,albsud,alb_bb,ith_bb,Tiso,Tabove
143	real ptoto,pcap,patm,airetot,ptotn,patmn,psea
144	! real ssum
145	character*1 yes
146	logical :: flagtset=.false. , flagps0=.false.
147	real val, val2, val3, val4 ! to store temporary variables
148
149	INTEGER :: itau
150
151	character(len=20) :: txt ! to store some text
152	character(len=50) :: surfacefile ! "surface.nc" (or equivalent file)
153	character(len=150) :: longline
154	integer :: count
155	real :: profile(llm+1) ! to store an atmospheric profile + surface value
156
157	! added by BC for equilibrium temperature startup
158	real teque
159
160	! added by RW for nuketharsis
161	real fact1
162	real fact2
163
164	! added by JVO for Titan
165	REAL tankCH4(ngridmx)
166
167	c sortie visu pour les champs dynamiques
168	c---------------------------------------
169	! INTEGER :: visuid
170	! real :: time_step,t_ops,t_wrt
171	! CHARACTER*80 :: visu_file
172
173	cpp = 0.
174	preff = 0.
175	pa = 0. ! to ensure disaster rather than minor error if we don`t
176	! make deliberate choice of these values elsewhere.
177
178	planet_type="titan"
179
180	! initialize "serial/parallel" related stuff
181	! (required because we call tabfi() below, before calling iniphysiq)
182	is_sequential=.true.
183	is_parallel=.false.
184	is_mpi_root=.true.
185	is_omp_root=.true.
186	is_master=.true.
187
188	! Load tracer number and names:
189	call infotrac_init
190	! allocate arrays
191	allocate(q(iip1,jjp1,llm,nqtot))
192	allocate(qsurf(ngridmx,nqtot))
193
194	! get value of nsoilmx and allocate corresponding arrays
195	! a default value of nsoilmx is set in comsoil_h
196	call getin_p("nsoilmx",nsoilmx)
197
198	allocate(tsoil(ngridmx,nsoilmx))
199	allocate(ith(iip1,jjp1,nsoilmx),ithfi(ngridmx,nsoilmx))
200
201	c=======================================================================
202	c Choice of the start file(s) to use
203	c=======================================================================
204	write(,) 'From which kind of files do you want to create new',
205	. 'start and startfi files'
206	write(,) ' 0 - from a file start_archive'
207	write(,) ' 1 - from files start and startfi'
208
209	c-----------------------------------------------------------------------
210	c Open file(s) to modify (start or start_archive)
211	c-----------------------------------------------------------------------
212
213	DO
214	read(,,iostat=ierr) choix_1
215	if ((choix_1 /= 0).OR.(choix_1 /=1)) EXIT
216	ENDDO
217
218	c Open start_archive
219	c ~~~~~~~~~~~~~~~~~~~~~~~~~~
220	if (choix_1.eq.0) then
221
222	write(,) 'Creating start files from:'
223	write(,) './start_archive.nc'
224	write(,)
225	fichnom = 'start_archive.nc'
226	ierr = NF_OPEN (fichnom, NF_NOWRITE,nid)
227	IF (ierr.NE.NF_NOERR) THEN
228	write(6,*)' Problem opening file:',fichnom
229	write(6,*)' ierr = ', ierr
230	CALL ABORT
231	ENDIF
232	tab0 = 50
233	Lmodif = 1
234
235	c OR open start and startfi files
236	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
237	else
238	write(,) 'Creating start files from:'
239	write(,) './start.nc and ./startfi.nc'
240	write(,)
241	fichnom = 'start.nc'
242	ierr = NF_OPEN (fichnom, NF_NOWRITE,nid_dyn)
243	IF (ierr.NE.NF_NOERR) THEN
244	write(6,*)' Problem opening file:',fichnom
245	write(6,*)' ierr = ', ierr
246	CALL ABORT
247	ENDIF
248
249	fichnom = 'startfi.nc'
250	ierr = NF_OPEN (fichnom, NF_NOWRITE,nid_fi)
251	IF (ierr.NE.NF_NOERR) THEN
252	write(6,*)' Problem opening file:',fichnom
253	write(6,*)' ierr = ', ierr
254	CALL ABORT
255	ENDIF
256
257	tab0 = 0
258	Lmodif = 0
259
260	endif
261
262
263	c=======================================================================
264	c INITIALISATIONS DIVERSES
265	c=======================================================================
266
267	! Initialize global tracer indexes (stored in tracer.h)
268	! ... this has to be done before phyetat0
269	call initracer(ngridmx,nqtot,tname)
270
271	! Take care of arrays in common modules
272	! ALLOCATE ARRAYS in surfdat_h (if not already done, e.g. when using start_archive)
273	IF (.not. ALLOCATED(albedodat)) ALLOCATE(albedodat(ngridmx))
274	IF (.not. ALLOCATED(phisfi)) ALLOCATE(phisfi(ngridmx))
275	IF (.not. ALLOCATED(zmea)) ALLOCATE(zmea(ngridmx))
276	IF (.not. ALLOCATED(zstd)) ALLOCATE(zstd(ngridmx))
277	IF (.not. ALLOCATED(zsig)) ALLOCATE(zsig(ngridmx))
278	IF (.not. ALLOCATED(zgam)) ALLOCATE(zgam(ngridmx))
279	IF (.not. ALLOCATED(zthe)) ALLOCATE(zthe(ngridmx))
280
281	c-----------------------------------------------------------------------
282	c Lecture du tableau des parametres du run (pour la dynamique)
283	c-----------------------------------------------------------------------
284	if (choix_1.eq.0) then
285
286	write(,) 'reading tab_cntrl START_ARCHIVE'
287	c
288	ierr = NF_INQ_VARID (nid, "controle", nvarid)
289	#ifdef NC_DOUBLE
290	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tab_cntrl)
291	#else
292	ierr = NF_GET_VAR_REAL(nid, nvarid, tab_cntrl)
293	#endif
294	c
295	else if (choix_1.eq.1) then
296
297	write(,) 'reading tab_cntrl START'
298	c
299	ierr = NF_INQ_VARID (nid_dyn, "controle", nvarid)
300	#ifdef NC_DOUBLE
301	ierr = NF_GET_VAR_DOUBLE(nid_dyn, nvarid, tab_cntrl)
302	#else
303	ierr = NF_GET_VAR_REAL(nid_dyn, nvarid, tab_cntrl)
304	#endif
305	c
306	write(,) 'reading tab_cntrl STARTFI'
307	c
308	ierr = NF_INQ_VARID (nid_fi, "controle", nvarid)
309	#ifdef NC_DOUBLE
310	ierr = NF_GET_VAR_DOUBLE(nid_fi, nvarid, tab_cntrl_bis)
311	#else
312	ierr = NF_GET_VAR_REAL(nid_fi, nvarid, tab_cntrl_bis)
313	#endif
314	c
315	do i=1,50
316	tab_cntrl(i+50)=tab_cntrl_bis(i)
317	enddo
318	write(,) 'printing tab_cntrl', tab_cntrl
319	do i=1,100
320	write(,) i,tab_cntrl(i)
321	enddo
322
323	write(,) 'Reading file START'
324	fichnom = 'start.nc'
325	CALL dynetat0(fichnom,vcov,ucov,teta,q,masse,
326	. ps,phis,time)
327
328	CALL iniconst
329	CALL inigeom
330
331	! Initialize the physics
332	CALL iniphysiq(iim,jjm,llm,
333	& (jjm-1)*iim+2,comm_lmdz,
334	& daysec,day_ini,dtphys,
335	& rlatu,rlatv,rlonu,rlonv,
336	& aire,cu,cv,rad,g,r,cpp,
337	& 1)
338
339	! Lmodif set to 0 to disable modifications possibility in phyeta0
340	Lmodif=0
341	write(,) 'Reading file STARTFI'
342	fichnom = 'startfi.nc'
343	CALL phyetat0(.true.,ngridmx,llm,fichnom,tab0,Lmodif,nsoilmx,
344	. nqtot,day_ini,time,
345	. tsurf,tsoil,emis,q2,qsurf,tankCH4)
346
347	! copy albedo and soil thermal inertia on (local) physics grid
348	do i=1,ngridmx
349	albfi(i) = albedodat(i)
350	do j=1,nsoilmx
351	ithfi(i,j) = inertiedat(i,j)
352	enddo
353	! build a surfithfi(:) using 1st layer of ithfi(:), which might
354	! be needed later on if reinitializing soil thermal inertia
355	surfithfi(i)=ithfi(i,1)
356	enddo
357	! also copy albedo and soil thermal inertia on (local) dynamics grid
358	! so that options below can manipulate either (but must then ensure
359	! to correctly recast things on physics grid)
360	call gr_fi_dyn(1,ngridmx,iip1,jjp1,albfi,alb)
361	call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
362	call gr_fi_dyn(1,ngridmx,iip1,jjp1,surfithfi,surfith)
363
364	endif
365	c-----------------------------------------------------------------------
366	c Initialisation des constantes dynamique
367	c-----------------------------------------------------------------------
368
369	kappa = tab_cntrl(9)
370	etot0 = tab_cntrl(12)
371	ptot0 = tab_cntrl(13)
372	ztot0 = tab_cntrl(14)
373	stot0 = tab_cntrl(15)
374	ang0 = tab_cntrl(16)
375	write(,) "Newstart: kappa,etot0,ptot0,ztot0,stot0,ang0"
376	write(,) kappa,etot0,ptot0,ztot0,stot0,ang0
377
378	! for vertical coordinate
379	preff=tab_cntrl(18) ! reference surface pressure
380	pa=tab_cntrl(17) ! reference pressure at which coord is purely pressure
381	!NB: in start_archive files tab_cntrl(17)=tab_cntrl(18)=0
382	write(,) "Newstart: preff=",preff," pa=",pa
383	yes=' '
384	do while ((yes.ne.'y').and.(yes.ne.'n'))
385	write(,) "Change the values of preff and pa? (y/n)"
386	read(*,fmt='(a)') yes
387	end do
388	if (yes.eq.'y') then
389	write(,)"New value of reference surface pressure preff?"
390	write(,)" (for Mars, typically preff=610)"
391	read(,) preff
392	write(,)"New value of reference pressure pa for purely"
393	write(,)"pressure levels (for hybrid coordinates)?"
394	write(,)" (for Mars, typically pa=20)"
395	read(,) pa
396	endif
397	c-----------------------------------------------------------------------
398	c Lecture du tab_cntrl et initialisation des constantes physiques
399	c - pour start: Lmodif = 0 => pas de modifications possibles
400	c (modif dans le tabfi de readfi + loin)
401	c - pour start_archive: Lmodif = 1 => modifications possibles
402	c-----------------------------------------------------------------------
403	if (choix_1.eq.0) then
404	! tabfi requires that input file be first opened by open_startphy(fichnom)
405	fichnom = 'start_archive.nc'
406	call open_startphy(fichnom)
407	call tabfi (ngridmx,nid,Lmodif,tab0,day_ini,lllm,p_rad,
408	. p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
409	else if (choix_1.eq.1) then
410	fichnom = 'startfi.nc'
411	call open_startphy(fichnom)
412	Lmodif=1 ! Lmodif set to 1 to allow modifications in phyeta0
413	call tabfi (ngridmx,nid_fi,Lmodif,tab0,day_ini,lllm,p_rad,
414	. p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
415	endif
416
417	if (p_omeg .eq. -9999.) then
418	p_omeg = 8.*atan(1.)/p_daysec
419	print*,"new value of omega",p_omeg
420	endif
421
422	rad = p_rad
423	omeg = p_omeg
424	g = p_g
425	cpp = p_cpp
426	mugaz = p_mugaz
427	daysec = p_daysec
428
429	if (p_omeg .eq. -9999.) then
430	p_omeg = 8.*atan(1.)/p_daysec
431	print*,"new value of omega",p_omeg
432	endif
433
434	kappa = 8.3141000./(p_mugazp_cpp) ! added by RDW
435
436	c=======================================================================
437	c INITIALISATIONS DIVERSES
438	c=======================================================================
439	! Load parameters from run.def file
440	CALL defrun_new( 99, .TRUE. )
441	! Initialize dynamics geometry and co. (which, when using start.nc may
442	! have changed because of modifications (tabi, preff, pa) above)
443	CALL iniconst
444	CALL inigeom
445	idum=-1
446	idum=0
447
448	! Initialize the physics for start_archive only
449	IF (choix_1.eq.0) THEN
450	CALL iniphysiq(iim,jjm,llm,
451	& (jjm-1)*iim+2,comm_lmdz,
452	& daysec,day_ini,dtphys,
453	& rlatu,rlatv,rlonu,rlonv,
454	& aire,cu,cv,rad,g,r,cpp,
455	& 1)
456	ENDIF
457
458	c=======================================================================
459	c lecture topographie, albedo, inertie thermique, relief sous-maille
460	c=======================================================================
461
462	if (choix_1.eq.0) then ! for start_archive files,
463	! where an external "surface.nc" file is needed
464
465	c do while((relief(1:3).ne.'mol').AND.(relief(1:3).ne.'pla'))
466	c write(,)
467	c write(,) 'choix du relief (mola,pla)'
468	c write(,) '(Topographie MGS MOLA, plat)'
469	c read(*,fmt='(a3)') relief
470	relief="mola"
471	c enddo
472
473	! First get the correct value of datadir, if not already done:
474	! default 'datadir' is set in "datafile_mod"
475	call getin_p("datadir",datadir)
476	write(,) 'Available surface data files are:'
477	filestring='ls '//trim(datadir)//'/'//
478	& trim(surfdir)//' \| grep .nc'
479	call system(filestring)
480	! but in ye old days these files were in datadir, so scan it as well
481	! for the sake of retro-compatibility
482	filestring='ls '//trim(datadir)//'/'//' \| grep .nc'
483	call system(filestring)
484
485	write(,) ''
486	write(,) 'Please choose the relevant file',
487	& ' (e.g. "surface_mars.nc")'
488	write(,) ' or "none" to not use any (and set planetary'
489	write(,) ' albedo and surface thermal inertia)'
490	read(*,fmt='(a50)') surfacefile
491
492	if (surfacefile.ne."none") then
493
494	CALL datareadnc(relief,surfacefile,phis,alb,surfith,
495	& zmeaS,zstdS,zsigS,zgamS,ztheS)
496	else
497	! specific case when not using a "surface.nc" file
498	phis(:,:)=0
499	zmeaS(:,:)=0
500	zstdS(:,:)=0
501	zsigS(:,:)=0
502	zgamS(:,:)=0
503	ztheS(:,:)=0
504
505	write(,) "Enter value of albedo of the bare ground:"
506	read(,) alb(1,1)
507	alb(:,:)=alb(1,1)
508
509	write(,) "Enter value of thermal inertia of soil:"
510	read(,) surfith(1,1)
511	surfith(:,:)=surfith(1,1)
512
513	endif ! of if (surfacefile.ne."none")
514
515	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
516	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,surfith,surfithfi)
517	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
518	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zmeaS,zmea)
519	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zstdS,zstd)
520	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zsigS,zsig)
521	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zgamS,zgam)
522	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ztheS,zthe)
523
524	endif ! of if (choix_1.eq.0)
525
526
527	c=======================================================================
528	c Lecture des fichiers (start ou start_archive)
529	c=======================================================================
530
531	if (choix_1.eq.0) then
532
533	write(,) 'Reading file START_ARCHIVE'
534	CALL lect_start_archive(ngridmx,llm,
535	& date,tsurf,tsoil,emis,q2,
536	& t,ucov,vcov,ps,teta,phisold_newgrid,q,qsurf,
537	& surfith,nid)
538	write(,) "OK, read start_archive file"
539	! copy soil thermal inertia
540	ithfi(:,:)=inertiedat(:,:)
541
542	ierr= NF_CLOSE(nid)
543
544	else if (choix_1.eq.1) then
545	!do nothing, start and startfi have already been read
546	else
547	CALL exit(1)
548	endif
549
550	dtvr = daysec/FLOAT(day_step)
551	dtphys = dtvr * FLOAT(iphysiq)
552
553	c=======================================================================
554	c
555	c=======================================================================
556
557	do ! infinite loop on list of changes
558
559	write(,)
560	write(,)
561	write(,) 'List of possible changes :'
562	write(,) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
563	write(,)
564	write(,) 'flat : no topography ("aquaplanet")'
565	write(,) 'set_ps_to_preff : used if changing preff with topo'
566	write(,) 'nuketharsis : no Tharsis bulge'
567	write(,) 'bilball : uniform albedo and thermal inertia'
568	write(,) 'coldspole : cold subsurface and high albedo at S.pole'
569	write(,) 'qname : change tracer name'
570	write(,) 't=profile : read temperature profile in profile.in'
571	write(,) 'q=0 : ALL tracer =zero'
572	write(,) 'q=x : give a specific uniform value to one tracer'
573	write(,) 'q=profile : specify a profile for a tracer'
574	! write(,) 'ini_q : tracers initialisation for chemistry, water an
575	! $d ice '
576	! write(,) 'ini_q-H2O : tracers initialisation for chemistry and
577	! $ice '
578	! write(,) 'ini_q-iceH2O : tracers initialisation for chemistry on
579	! $ly '
580	write(,) 'noglacier : Remove tropical H2O ice if \|lat\|<45'
581	write(,) 'watercapn : H20 ice on permanent N polar cap '
582	write(,) 'watercaps : H20 ice on permanent S polar cap '
583	write(,) 'noacglac : H2O ice across Noachis Terra'
584	write(,) 'oborealis : H2O ice across Vastitas Borealis'
585	write(,) 'iceball : Thick ice layer all over surface'
586	write(,) 'supercontinent: Create a continent of given Ab and TI'
587	write(,) 'wetstart : start with a wet atmosphere'
588	write(,) 'isotherm : Isothermal Temperatures, wind set to zero'
589	write(,) 'radequi : Earth-like radiative equilibrium temperature
590	$ profile (lat-alt) and winds set to zero'
591	write(,) 'coldstart : Start X K above the CO2 frost point and
592	$set wind to zero (assumes 100% CO2)'
593	write(,) 'co2ice=0 : remove CO2 polar cap'
594	write(,) 'ptot : change total pressure'
595	write(,) 'emis : change surface emissivity'
596	write(,) 'therm_ini_s : Set soil thermal inertia to reference sur
597	&face values'
598	write(,) 'slab_ocean_0 : initialisation of slab
599	$ocean variables'
600
601	write(,)
602	write(,) 'Change to perform ?'
603	write(,) ' (enter keyword or return to end)'
604	write(,)
605
606	read(*,fmt='(a20)') modif
607	if (modif(1:1) .eq. ' ')then
608	exit ! exit loop on changes
609	endif
610
611	write(,)
612	write(,) trim(modif) , ' : '
613
614	c 'flat : no topography ("aquaplanet")'
615	c -------------------------------------
616	if (trim(modif) .eq. 'flat') then
617	c set topo to zero
618	z_reel(1:iip1,1:jjp1)=0
619	phis(1:iip1,1:jjp1)=g*z_reel(1:iip1,1:jjp1)
620	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
621	write(,) 'topography set to zero.'
622	write(,) 'WARNING : the subgrid topography parameters',
623	& ' were not set to zero ! => set calllott to F'
624
625	c Choice of surface pressure
626	yes=' '
627	do while ((yes.ne.'y').and.(yes.ne.'n'))
628	write(,) 'Do you wish to choose homogeneous surface',
629	& 'pressure (y) or let newstart interpolate ',
630	& ' the previous field (n)?'
631	read(*,fmt='(a)') yes
632	end do
633	if (yes.eq.'y') then
634	flagps0=.true.
635	write(,) 'New value for ps (Pa) ?'
636	201 read(,,iostat=ierr) patm
637	if(ierr.ne.0) goto 201
638	write(,) patm
639	if (patm.eq.-9999.) then
640	patm = preff
641	write(,) "we set patm = preff", preff
642	endif
643	write(,)
644	write(,) ' new ps everywhere (Pa) = ', patm
645	write(,)
646	do j=1,jjp1
647	do i=1,iip1
648	ps(i,j)=patm
649	enddo
650	enddo
651	end if
652
653	c 'set_ps_to_preff' : used if changing preff with topo
654	c ----------------------------------------------------
655	else if (trim(modif) .eq. 'set_ps_to_preff') then
656	do j=1,jjp1
657	do i=1,iip1
658	ps(i,j)=preff
659	enddo
660	enddo
661
662	c 'nuketharsis : no tharsis bulge for Early Mars'
663	c ---------------------------------------------
664	else if (trim(modif) .eq. 'nuketharsis') then
665
666	DO j=1,jjp1
667	DO i=1,iim
668	ig=1+(j-2)*iim +i
669	if(j.eq.1) ig=1
670	if(j.eq.jjp1) ig=ngridmx
671
672	fact1=(((rlonv(i)180./pi)+100)*2 +
673	& (rlatu(j)180./pi)2)/65*2
674	fact2=exp( -fact1**2.5 )
675
676	phis(i,j) = phis(i,j) - (phis(i,j)+4000.g)fact2
677
678	! if(phis(i,j).gt.2500.*g)then
679	! if(rlatu(j)*180./pi.gt.-80.)then ! avoid chopping south polar cap
680	! phis(i,j)=2500.*g
681	! endif
682	! endif
683
684	ENDDO
685	ENDDO
686	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
687
688
689	c bilball : uniform albedo, thermal inertia
690	c -----------------------------------------
691	else if (trim(modif) .eq. 'bilball') then
692	write(,) 'constante albedo and iner.therm:'
693	write(,) 'New value for albedo (ex: 0.25) ?'
694	101 read(,,iostat=ierr) alb_bb
695	if(ierr.ne.0) goto 101
696	write(,)
697	write(,) ' uniform albedo (new value):',alb_bb
698	write(,)
699
700	write(,) 'New value for thermal inertia (eg: 247) ?'
701	102 read(,,iostat=ierr) ith_bb
702	if(ierr.ne.0) goto 102
703	write(,) 'uniform thermal inertia (new value):',ith_bb
704	DO j=1,jjp1
705	DO i=1,iip1
706	alb(i,j) = alb_bb ! albedo
707	do isoil=1,nsoilmx
708	ith(i,j,isoil) = ith_bb ! thermal inertia
709	enddo
710	END DO
711	END DO
712	! CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ith,ithfi)
713	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
714	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
715
716	c coldspole : sous-sol de la calotte sud toujours froid
717	c -----------------------------------------------------
718	else if (trim(modif) .eq. 'coldspole') then
719	write(,)'new value for the subsurface temperature',
720	& ' beneath the permanent southern polar cap ? (eg: 141 K)'
721	103 read(,,iostat=ierr) tsud
722	if(ierr.ne.0) goto 103
723	write(,)
724	write(,) ' new value of the subsurface temperature:',tsud
725	c nouvelle temperature sous la calotte permanente
726	do l=2,nsoilmx
727	tsoil(ngridmx,l) = tsud
728	end do
729
730
731	write(,)'new value for the albedo',
732	& 'of the permanent southern polar cap ? (eg: 0.75)'
733	104 read(,,iostat=ierr) albsud
734	if(ierr.ne.0) goto 104
735	write(,)
736
737	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
738	c Option 1: only the albedo of the pole is modified :
739	albfi(ngridmx)=albsud
740	write(,) 'ig=',ngridmx,' albedo perennial cap ',
741	& albfi(ngridmx)
742
743	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
744	c Option 2 A haute resolution : coordonnee de la vrai calotte ~
745	c DO j=1,jjp1
746	c DO i=1,iip1
747	c ig=1+(j-2)*iim +i
748	c if(j.eq.1) ig=1
749	c if(j.eq.jjp1) ig=ngridmx
750	c if ((rlatu(j)*180./pi.lt.-84.).and.
751	c & (rlatu(j)*180./pi.gt.-91.).and.
752	c & (rlonv(i)*180./pi.gt.-91.).and.
753	c & (rlonv(i)*180./pi.lt.0.)) then
754	cc albedo de la calotte permanente fixe a albsud
755	c alb(i,j)=albsud
756	c write(,) 'lat=',rlatu(j)*180./pi,
757	c & ' lon=',rlonv(i)*180./pi
758	cc fin de la condition sur les limites de la calotte permanente
759	c end if
760	c ENDDO
761	c ENDDO
762	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
763
764	c CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
765
766
767	c ptot : Modification of the total pressure: ice + current atmosphere
768	c -------------------------------------------------------------------
769	else if (trim(modif).eq.'ptot') then
770
771	c calcul de la pression totale glace + atm actuelle
772	patm=0.
773	airetot=0.
774	pcap=0.
775	DO j=1,jjp1
776	DO i=1,iim
777	ig=1+(j-2)*iim +i
778	if(j.eq.1) ig=1
779	if(j.eq.jjp1) ig=ngridmx
780	patm = patm + ps(i,j)*aire(i,j)
781	airetot= airetot + aire(i,j)
782	ENDDO
783	ENDDO
784	ptoto = pcap + patm
785
786	print*,'Current total pressure at surface (atm) ',
787	& ptoto/airetot
788
789	print*,'new value?'
790	read(,) ptotn
791	ptotn=ptotn*airetot
792	patmn=ptotn-pcap
793	print*,'ptoto,patm,ptotn,patmn'
794	print*,ptoto,patm,ptotn,patmn
795	print*,'Mult. factor for pressure (atm only)', patmn/patm
796	do j=1,jjp1
797	do i=1,iip1
798	ps(i,j)=ps(i,j)*patmn/patm
799	enddo
800	enddo
801
802
803
804	c Correction pour la conservation des traceurs
805	yes=' '
806	do while ((yes.ne.'y').and.(yes.ne.'n'))
807	write(,) 'Do you wish to conserve tracer total mass (y)',
808	& ' or tracer mixing ratio (n) ?'
809	read(*,fmt='(a)') yes
810	end do
811
812	if (yes.eq.'y') then
813	write(,) 'OK : conservation of tracer total mass'
814	DO iq =1, nqtot
815	DO l=1,llm
816	DO j=1,jjp1
817	DO i=1,iip1
818	q(i,j,l,iq)=q(i,j,l,iq)*patm/patmn
819	ENDDO
820	ENDDO
821	ENDDO
822	ENDDO
823	else
824	write(,) 'OK : conservation of tracer mixing ratio'
825	end if
826
827	c Correction pour la pression au niveau de la mer
828	yes=' '
829	do while ((yes.ne.'y').and.(yes.ne.'n'))
830	write(,) 'Do you wish fix pressure at sea level (y)',
831	& ' or not (n) ?'
832	read(*,fmt='(a)') yes
833	end do
834
835	if (yes.eq.'y') then
836	write(,) 'Value?'
837	read(,,iostat=ierr) psea
838	DO i=1,iip1
839	DO j=1,jjp1
840	ps(i,j)=psea
841
842	ENDDO
843	ENDDO
844	write(,) 'psea=',psea
845	else
846	write(,) 'no'
847	end if
848
849
850	c emis : change surface emissivity (added by RW)
851	c ----------------------------------------------
852	else if (trim(modif).eq.'emis') then
853
854	print*,'new value?'
855	read(,) emisread
856
857	do i=1,ngridmx
858	emis(i)=emisread
859	enddo
860
861	c qname : change tracer name
862	c --------------------------
863	else if (trim(modif).eq.'qname') then
864	yes='y'
865	do while (yes.eq.'y')
866	write(,) 'Which tracer name do you want to change ?'
867	do iq=1,nqtot
868	write(*,'(i3,a3,a20)')iq,' : ',trim(tname(iq))
869	enddo
870	write(*,'(a35,i3)')
871	& '(enter tracer number; between 1 and ',nqtot
872	write(,)' or any other value to quit this option)'
873	read(,) iq
874	if ((iq.ge.1).and.(iq.le.nqtot)) then
875	write(,)'Change tracer name ',trim(tname(iq)),' to ?'
876	read(,) txt
877	tname(iq)=txt
878	write(,)'Do you want to change another tracer name (y/n)?'
879	read(*,'(a)') yes
880	else
881	! inapropiate value of iq; quit this option
882	yes='n'
883	endif ! of if ((iq.ge.1).and.(iq.le.nqtot))
884	enddo ! of do while (yes.ne.'y')
885
886	c q=0 : set tracers to zero
887	c -------------------------
888	else if (trim(modif).eq.'q=0') then
889	c mise a 0 des q (traceurs)
890	write(,) 'Tracers set to 0 (1.E-30 in fact)'
891	DO iq =1, nqtot
892	DO l=1,llm
893	DO j=1,jjp1
894	DO i=1,iip1
895	q(i,j,l,iq)=1.e-30
896	ENDDO
897	ENDDO
898	ENDDO
899	ENDDO
900
901	c set surface tracers to zero
902	DO iq =1, nqtot
903	DO ig=1,ngridmx
904	qsurf(ig,iq)=0.
905	ENDDO
906	ENDDO
907
908	c q=x : initialise tracer manually
909	c --------------------------------
910	else if (trim(modif).eq.'q=x') then
911	write(,) 'Which tracer do you want to modify ?'
912	do iq=1,nqtot
913	write(,)iq,' : ',trim(tname(iq))
914	enddo
915	write(,) '(choose between 1 and ',nqtot,')'
916	read(,) iq
917	write(,)'mixing ratio of tracer ',trim(tname(iq)),
918	& ' ? (kg/kg)'
919	read(,) val
920	DO l=1,llm
921	DO j=1,jjp1
922	DO i=1,iip1
923	q(i,j,l,iq)=val
924	ENDDO
925	ENDDO
926	ENDDO
927	write(,) 'SURFACE value of tracer ',trim(tname(iq)),
928	& ' ? (kg/m2)'
929	read(,) val
930	DO ig=1,ngridmx
931	qsurf(ig,iq)=val
932	ENDDO
933
934	c t=profile : initialize temperature with a given profile
935	c -------------------------------------------------------
936	else if (trim(modif) .eq. 't=profile') then
937	write(,) 'Temperature profile from ASCII file'
938	write(,) "'profile.in' e.g. 1D output"
939	write(,) "(one value per line in file; starting with"
940	write(,) "surface value, the 1st atmospheric layer"
941	write(,) "followed by 2nd, etc. up to top of atmosphere)"
942	txt="profile.in"
943	open(33,file=trim(txt),status='old',form='formatted',
944	& iostat=ierr)
945	if (ierr.eq.0) then
946	! OK, found file 'profile_...', load the profile
947	do l=1,llm+1
948	read(33,*,iostat=ierr) profile(l)
949	write(,) profile(l)
950	if (ierr.ne.0) then ! something went wrong
951	exit ! quit loop
952	endif
953	enddo
954	if (ierr.eq.0) then
955	tsurf(1:ngridmx)=profile(1)
956	tsoil(1:ngridmx,1:nsoilmx)=profile(1)
957	do l=1,llm
958	Tset(1:iip1,1:jjp1,l)=profile(l+1)
959	flagtset=.true.
960	enddo
961	ucov(1:iip1,1:jjp1,1:llm)=0.
962	vcov(1:iip1,1:jjm,1:llm)=0.
963	q2(1:ngridmx,1:llm+1)=0.
964	else
965	write(,)'problem reading file ',trim(txt),' !'
966	write(,)'No modifications to temperature'
967	endif
968	else
969	write(,)'Could not find file ',trim(txt),' !'
970	write(,)'No modifications to temperature'
971	endif
972
973	c q=profile : initialize tracer with a given profile
974	c --------------------------------------------------
975	else if (trim(modif) .eq. 'q=profile') then
976	write(,) 'Tracer profile will be sought in ASCII file'
977	write(,) "'profile_tracer' where 'tracer' is tracer name"
978	write(,) "(one value per line in file; starting with"
979	write(,) "surface value, the 1st atmospheric layer"
980	write(,) "followed by 2nd, etc. up to top of atmosphere)"
981	write(,) 'Which tracer do you want to set?'
982	do iq=1,nqtot
983	write(,)iq,' : ',trim(tname(iq))
984	enddo
985	write(,) '(choose between 1 and ',nqtot,')'
986	read(,) iq
987	if ((iq.lt.1).or.(iq.gt.nqtot)) then
988	! wrong value for iq, go back to menu
989	write(,) "wrong input value:",iq
990	cycle
991	endif
992	! look for input file 'profile_tracer'
993	txt="profile_"//trim(tname(iq))
994	open(41,file=trim(txt),status='old',form='formatted',
995	& iostat=ierr)
996	if (ierr.eq.0) then
997	! OK, found file 'profile_...', load the profile
998	do l=1,llm+1
999	read(41,*,iostat=ierr) profile(l)
1000	if (ierr.ne.0) then ! something went wrong
1001	exit ! quit loop
1002	endif
1003	enddo
1004	if (ierr.eq.0) then
1005	! initialize tracer values
1006	qsurf(:,iq)=profile(1)
1007	do l=1,llm
1008	q(:,:,l,iq)=profile(l+1)
1009	enddo
1010	write(,)'OK, tracer ',trim(tname(iq)),' initialized '
1011	& ,'using values from file ',trim(txt)
1012	else
1013	write(,)'problem reading file ',trim(txt),' !'
1014	write(,)'No modifications to tracer ',trim(tname(iq))
1015	endif
1016	else
1017	write(,)'Could not find file ',trim(txt),' !'
1018	write(,)'No modifications to tracer ',trim(tname(iq))
1019	endif
1020
1021
1022	c wetstart : wet atmosphere with a north to south gradient
1023	c --------------------------------------------------------
1024	else if (trim(modif) .eq. 'wetstart') then
1025	! check that there is indeed a water vapor tracer
1026
1027	write(,) "No water vapour tracer! Can't use this option"
1028	stop
1029
1030	c noglacier : remove tropical water ice (to initialize high res sim)
1031	c --------------------------------------------------
1032	else if (trim(modif) .eq. 'noglacier') then
1033
1034	write(,) "No water ice tracer! Can't use this option"
1035	stop
1036
1037
1038
1039	c watercapn : H20 ice on permanent northern cap
1040	c --------------------------------------------------
1041	else if (trim(modif) .eq. 'watercapn') then
1042
1043	write(,) "No water ice tracer! Can't use this option"
1044	stop
1045
1046	c watercaps : H20 ice on permanent southern cap
1047	c -------------------------------------------------
1048	else if (trim(modif) .eq. 'watercaps') then
1049
1050	write(,) "No water ice tracer! Can't use this option"
1051	stop
1052
1053	c noacglac : H2O ice across highest terrain
1054	c --------------------------------------------
1055	else if (trim(modif) .eq. 'noacglac') then
1056
1057	write(,) "No water ice tracer! Can't use this option"
1058	stop
1059
1060	c oborealis : H2O oceans across Vastitas Borealis
1061	c -----------------------------------------------
1062	else if (trim(modif) .eq. 'oborealis') then
1063
1064	write(,) "No water ice tracer! Can't use this option"
1065	stop
1066
1067	c iborealis : H2O ice in Northern plains
1068	c --------------------------------------
1069	else if (trim(modif) .eq. 'iborealis') then
1070
1071	write(,) "No water ice tracer! Can't use this option"
1072	stop
1073
1074	c oceanball : H2O liquid everywhere
1075	c ----------------------------
1076	else if (trim(modif) .eq. 'oceanball') then
1077
1078	write(,) "No water ice tracer! Can't use this option"
1079	stop
1080
1081	c iceball : H2O ice everywhere
1082	c ----------------------------
1083	else if (trim(modif) .eq. 'iceball') then
1084
1085	write(,) "No water ice tracer! Can't use this option"
1086	stop
1087
1088	c supercontinent : H2O ice everywhere
1089	c ----------------------------
1090	else if (trim(modif) .eq. 'supercontinent') then
1091
1092	write(,) "No water ice tracer! Can't use this option"
1093	stop
1094
1095	c isotherm : Isothermal temperatures and no winds
1096	c -----------------------------------------------
1097	else if (trim(modif) .eq. 'isotherm') then
1098
1099	write(,)'Isothermal temperature of the atmosphere,
1100	& surface and subsurface'
1101	write(,) 'Value of this temperature ? :'
1102	203 read(,,iostat=ierr) Tiso
1103	if(ierr.ne.0) goto 203
1104
1105	tsurf(1:ngridmx)=Tiso
1106
1107	tsoil(1:ngridmx,1:nsoilmx)=Tiso
1108
1109	Tset(1:iip1,1:jjp1,1:llm)=Tiso
1110	flagtset=.true.
1111
1112	t(1:iip1,1:jjp1,1:llm)=Tiso
1113	!! otherwise hydrost. integrations below
1114	!! use the wrong temperature
1115	!! -- NB: Tset might be useless
1116
1117	ucov(1:iip1,1:jjp1,1:llm)=0
1118	vcov(1:iip1,1:jjm,1:llm)=0
1119	q2(1:ngridmx,1:llm+1)=0
1120
1121	c radequi : Radiative equilibrium profile of temperatures and no winds
1122	c --------------------------------------------------------------------
1123	else if (trim(modif) .eq. 'radequi') then
1124
1125	write(,)'radiative equilibrium temperature profile'
1126
1127	do ig=1, ngridmx
1128	teque= 335.0-60.0sin(latfi(ig))sin(latfi(ig))-
1129	& 10.0cos(latfi(ig))cos(latfi(ig))
1130	tsurf(ig) = MAX(220.0,teque)
1131	end do
1132	do l=2,nsoilmx
1133	do ig=1, ngridmx
1134	tsoil(ig,l) = tsurf(ig)
1135	end do
1136	end do
1137	DO j=1,jjp1
1138	DO i=1,iim
1139	Do l=1,llm
1140	teque=335.-60.0sin(rlatu(j))sin(rlatu(j))-
1141	& 10.0cos(rlatu(j))cos(rlatu(j))
1142	Tset(i,j,l)=MAX(220.0,teque)
1143	end do
1144	end do
1145	end do
1146	flagtset=.true.
1147	ucov(1:iip1,1:jjp1,1:llm)=0
1148	vcov(1:iip1,1:jjm,1:llm)=0
1149	q2(1:ngridmx,1:llm+1)=0
1150
1151	c coldstart : T set 1K above CO2 frost point and no winds
1152	c ------------------------------------------------
1153	else if (trim(modif) .eq. 'coldstart') then
1154
1155	write(,)'set temperature of the atmosphere,'
1156	&,'surface and subsurface how many degrees above CO2 frost point?'
1157	204 read(,,iostat=ierr) Tabove
1158	if(ierr.ne.0) goto 204
1159
1160	DO j=1,jjp1
1161	DO i=1,iim
1162	ig=1+(j-2)*iim +i
1163	if(j.eq.1) ig=1
1164	if(j.eq.jjp1) ig=ngridmx
1165	tsurf(ig) = (-3167.8)/(log(.01*ps(i,j))-23.23)+Tabove
1166	END DO
1167	END DO
1168	do l=1,nsoilmx
1169	do ig=1, ngridmx
1170	tsoil(ig,l) = tsurf(ig)
1171	end do
1172	end do
1173	DO j=1,jjp1
1174	DO i=1,iim
1175	Do l=1,llm
1176	pp = aps(l) +bps(l)*ps(i,j)
1177	Tset(i,j,l)=(-3167.8)/(log(.01*pp)-23.23)+Tabove
1178	end do
1179	end do
1180	end do
1181
1182	flagtset=.true.
1183	ucov(1:iip1,1:jjp1,1:llm)=0
1184	vcov(1:iip1,1:jjm,1:llm)=0
1185	q2(1:ngridmx,1:llm+1)=0
1186
1187
1188	c co2ice=0 : remove CO2 polar ice caps'
1189	c ------------------------------------------------
1190	else if (trim(modif) .eq. 'co2ice=0') then
1191	write(,) "Can't remove CO2 ice!! (no co2_ice tracer)"
1192
1193
1194	! therm_ini_s: (re)-set soil thermal inertia to reference surface values
1195	! ----------------------------------------------------------------------
1196
1197	else if (trim(modif) .eq. 'therm_ini_s') then
1198	! write(,)"surfithfi(1):",surfithfi(1)
1199	do isoil=1,nsoilmx
1200	inertiedat(1:ngridmx,isoil)=surfithfi(1:ngridmx)
1201	enddo
1202	write(,)'OK: Soil thermal inertia has been reset to referenc
1203	&e surface values'
1204	! write(,)"inertiedat(1,1):",inertiedat(1,1)
1205	ithfi(:,:)=inertiedat(:,:)
1206	! recast ithfi() onto ith()
1207	call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
1208	! Check:
1209	! do i=1,iip1
1210	! do j=1,jjp1
1211	! do isoil=1,nsoilmx
1212	! write(77,*) i,j,isoil," ",ith(i,j,isoil)
1213	! enddo
1214	! enddo
1215	! enddo
1216
1217
1218
1219	c slab_ocean_initialisation
1220	c ------------------------------------------------
1221	else if (trim(modif) .eq. 'slab_ocean_0') then
1222
1223	write(,) "No ocean yet on Titan ! Can't use this option"
1224	stop
1225
1226	else
1227	write(,) ' Unknown (misspelled?) option!!!'
1228	end if ! of if (trim(modif) .eq. '...') elseif ...
1229
1230
1231
1232	enddo ! of do ! infinite loop on liste of changes
1233
1234	999 continue
1235
1236
1237	c================================================================
1238	c Initialisation for methane surface tank (added by JVO 2017)
1239	c================================================================
1240	DO ig=1, ngridmx
1241	tankCH4(ig)=2.0
1242	ENDDO
1243
1244	c=======================================================================
1245	c Correct pressure on the new grid (menu 0)
1246	c=======================================================================
1247
1248
1249	if ((choix_1.eq.0).and.(.not.flagps0)) then
1250	r = 1000.*8.31/mugaz
1251
1252	do j=1,jjp1
1253	do i=1,iip1
1254	ps(i,j) = ps(i,j) *
1255	. exp((phisold_newgrid(i,j)-phis(i,j)) /
1256	. (t(i,j,1) * r))
1257	end do
1258	end do
1259
1260	c periodicite de ps en longitude
1261	do j=1,jjp1
1262	ps(1,j) = ps(iip1,j)
1263	end do
1264	end if
1265
1266
1267	c=======================================================================
1268	c=======================================================================
1269
1270	c=======================================================================
1271	c Initialisation de la physique / ecriture de newstartfi :
1272	c=======================================================================
1273
1274
1275	CALL inifilr
1276	CALL pression(ip1jmp1, ap, bp, ps, p3d)
1277
1278	c-----------------------------------------------------------------------
1279	c Initialisation pks:
1280	c-----------------------------------------------------------------------
1281
1282	CALL exner_hyb(ip1jmp1, ps, p3d, beta, pks, pk, pkf)
1283	! Calcul de la temperature potentielle teta
1284
1285	if (flagtset) then
1286	DO l=1,llm
1287	DO j=1,jjp1
1288	DO i=1,iim
1289	teta(i,j,l) = Tset(i,j,l) * cpp/pk(i,j,l)
1290	ENDDO
1291	teta (iip1,j,l)= teta (1,j,l)
1292	ENDDO
1293	ENDDO
1294	else if (choix_1.eq.0) then
1295	DO l=1,llm
1296	DO j=1,jjp1
1297	DO i=1,iim
1298	teta(i,j,l) = t(i,j,l) * cpp/pk(i,j,l)
1299	ENDDO
1300	teta (iip1,j,l)= teta (1,j,l)
1301	ENDDO
1302	ENDDO
1303	endif
1304
1305	C Calcul intermediaire
1306	c
1307	if (choix_1.eq.0) then
1308	CALL massdair( p3d, masse )
1309	c
1310	print *,' ALPHAX ',alphax
1311
1312	DO l = 1, llm
1313	DO i = 1, iim
1314	xppn(i) = aire( i, 1 ) * masse( i , 1 , l )
1315	xpps(i) = aire( i,jjp1 ) * masse( i , jjp1 , l )
1316	ENDDO
1317	xpn = SUM(xppn)/apoln
1318	xps = SUM(xpps)/apols
1319	DO i = 1, iip1
1320	masse( i , 1 , l ) = xpn
1321	masse( i , jjp1 , l ) = xps
1322	ENDDO
1323	ENDDO
1324	endif
1325	phis(iip1,:) = phis(1,:)
1326
1327	itau=0
1328	if (choix_1.eq.0) then
1329	day_ini=int(date)
1330	endif
1331	c
1332	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
1333
1334	CALL caldyn0( itau,ucov,vcov,teta,ps,masse,pk,phis ,
1335	* phi,w, pbaru,pbarv,day_ini+time )
1336
1337
1338	CALL dynredem0("restart.nc",day_ini,phis)
1339	CALL dynredem1("restart.nc",0.0,vcov,ucov,teta,q,masse,ps)
1340	C
1341	C Ecriture etat initial physique
1342	C
1343
1344	call physdem0("restartfi.nc",lonfi,latfi,nsoilmx,ngridmx,llm,
1345	& nqtot,dtphys,real(day_ini),0.0,
1346	& airefi,albfi,ithfi,zmea,zstd,zsig,zgam,zthe)
1347	call physdem1("restartfi.nc",nsoilmx,ngridmx,llm,nqtot,
1348	& dtphys,real(day_ini),
1349	& tsurf,tsoil,emis,q2,qsurf,tankCH4)
1350
1351	c=======================================================================
1352	c Formats
1353	c=======================================================================
1354
1355	1 FORMAT(//10x,'la valeur de im =',i4,2x,'lue sur le fichier de dema
1356	*rrage est differente de la valeur parametree iim =',i4//)
1357	2 FORMAT(//10x,'la valeur de jm =',i4,2x,'lue sur le fichier de dema
1358	*rrage est differente de la valeur parametree jjm =',i4//)
1359	3 FORMAT(//10x,'la valeur de lllm =',i4,2x,'lue sur le fichier demar
1360	*rage est differente de la valeur parametree llm =',i4//)
1361
1362	write(,) "newstart: All is well that ends well."
1363
1364	end
1365

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

Download in other formats:

Original Format