Context Navigation

newstart.F @ 1647

Last change on this file since 1647 was 1647, checked in by jvatant, 8 years ago

+ Major clean of the new LMDZ.TITAN from too-generic options and routines (water, co2, ocean, surface type ...)
+ From this revision LMDZ.TITAN begins to be really separated from LMDZ.GENERIC
+ Partial desactivation of aerosols, only the dummy case is still enabled to keep the code running ( new aerosol routines to come in followings commits )

JVO

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

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

Download in other formats:

Original Format