Context Navigation

newstart.F @ 1889

Last change on this file since 1889 was 1889, checked in by jvatant, 7 years ago
Making chemistry handling more flexible - Step 2 + Added the calculation of the pressure grid in newstart using Vervack profile in gr_kim_vervack routine + Next step : regridding ! --JVO
File size: 40.5 KB

Line
1	C======================================================================
2	PROGRAM newstart
3	c=======================================================================
4	c
5	c
6	c Auteur: Christophe Hourdin/Francois Forget/Yann Wanherdrick
7	c ------
8	c
9	c Objet: Create or modify the initial state for the LMD Titan GCM
10	c ----- (fichiers NetCDF start et startfi)
11	c
12	c
13	c=======================================================================
14
15	use mod_phys_lmdz_para, only: is_parallel, is_sequential,
16	& is_mpi_root, is_omp_root,
17	& is_master
18	use infotrac, only: infotrac_init, nqtot, tname
19	USE comchem_h, ONLY: nlaykim_up, nlaykimold
20	USE comsoil_h, ONLY: nsoilmx, layer, mlayer, inertiedat
21	USE surfdat_h, ONLY: phisfi, albedodat,
22	& zmea, zstd, zsig, zgam, zthe
23	use datafile_mod, only: datadir, surfdir
24	use ioipsl_getin_p_mod, only: getin_p
25	use control_mod, only: day_step, iphysiq, anneeref, planet_type
26	use phyredem, only: physdem0, physdem1
27	use iostart, only: open_startphy
28	use filtreg_mod, only: inifilr
29	USE mod_const_mpi, ONLY: COMM_LMDZ
30	USE comvert_mod, ONLY: ap,bp,aps,bps,pa,preff
31	USE comconst_mod, ONLY: lllm,daysec,dtvr,dtphys,cpp,kappa,
32	. rad,omeg,g,r,pi
33	USE serre_mod, ONLY: alphax
34	USE temps_mod, ONLY: day_ini
35	USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
36	use tabfi_mod, only: tabfi
37	use iniphysiq_mod, only: iniphysiq
38	use phyetat0_mod, only: phyetat0
39	use tracer_h
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	c=======================================================================
263	c INITIALISATIONS DIVERSES
264	c=======================================================================
265
266	! Initialize global tracer indexes (stored in tracer.h)
267	! ... this has to be done before phyetat0
268	call initracer2(nqtot,tname)
269
270	! Take care of arrays in common modules
271	! ALLOCATE ARRAYS in surfdat_h (if not already done, e.g. when using start_archive)
272	IF (.not. ALLOCATED(albedodat)) ALLOCATE(albedodat(ngridmx))
273	IF (.not. ALLOCATED(phisfi)) ALLOCATE(phisfi(ngridmx))
274	IF (.not. ALLOCATED(zmea)) ALLOCATE(zmea(ngridmx))
275	IF (.not. ALLOCATED(zstd)) ALLOCATE(zstd(ngridmx))
276	IF (.not. ALLOCATED(zsig)) ALLOCATE(zsig(ngridmx))
277	IF (.not. ALLOCATED(zgam)) ALLOCATE(zgam(ngridmx))
278	IF (.not. ALLOCATED(zthe)) ALLOCATE(zthe(ngridmx))
279
280	c-----------------------------------------------------------------------
281	c Lecture du tableau des parametres du run (pour la dynamique)
282	c-----------------------------------------------------------------------
283	if (choix_1.eq.0) then
284
285	write(,) 'reading tab_cntrl START_ARCHIVE'
286	c
287	ierr = NF_INQ_VARID (nid, "controle", nvarid)
288	#ifdef NC_DOUBLE
289	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tab_cntrl)
290	#else
291	ierr = NF_GET_VAR_REAL(nid, nvarid, tab_cntrl)
292	#endif
293	c
294	else if (choix_1.eq.1) then
295
296	write(,) 'reading tab_cntrl START'
297	c
298	ierr = NF_INQ_VARID (nid_dyn, "controle", nvarid)
299	#ifdef NC_DOUBLE
300	ierr = NF_GET_VAR_DOUBLE(nid_dyn, nvarid, tab_cntrl)
301	#else
302	ierr = NF_GET_VAR_REAL(nid_dyn, nvarid, tab_cntrl)
303	#endif
304	c
305	write(,) 'reading tab_cntrl STARTFI'
306	c
307	ierr = NF_INQ_VARID (nid_fi, "controle", nvarid)
308	#ifdef NC_DOUBLE
309	ierr = NF_GET_VAR_DOUBLE(nid_fi, nvarid, tab_cntrl_bis)
310	#else
311	ierr = NF_GET_VAR_REAL(nid_fi, nvarid, tab_cntrl_bis)
312	#endif
313	c
314	do i=1,50
315	tab_cntrl(i+50)=tab_cntrl_bis(i)
316	enddo
317	write(,) 'printing tab_cntrl', tab_cntrl
318	do i=1,100
319	write(,) i,tab_cntrl(i)
320	enddo
321
322	write(,) 'Reading file START'
323	fichnom = 'start.nc'
324	CALL dynetat0(fichnom,vcov,ucov,teta,q,masse,
325	. ps,phis,time)
326
327	CALL iniconst
328	CALL inigeom
329
330	! Initialize the physics
331	CALL iniphysiq(iim,jjm,llm,
332	& (jjm-1)*iim+2,comm_lmdz,
333	& daysec,day_ini,dtphys,
334	& rlatu,rlatv,rlonu,rlonv,
335	& aire,cu,cv,rad,g,r,cpp,
336	& 1)
337
338	! Lmodif set to 0 to disable modifications possibility in phyeta0
339	Lmodif=0
340	write(,) 'Reading file STARTFI'
341	fichnom = 'startfi.nc'
342	CALL phyetat0(.true.,ngridmx,llm,fichnom,tab0,Lmodif,nsoilmx,
343	. nqtot,day_ini,time,
344	. tsurf,tsoil,emis,q2,qsurf,tankCH4)
345
346	! copy albedo and soil thermal inertia on (local) physics grid
347	do i=1,ngridmx
348	albfi(i) = albedodat(i)
349	do j=1,nsoilmx
350	ithfi(i,j) = inertiedat(i,j)
351	enddo
352	! build a surfithfi(:) using 1st layer of ithfi(:), which might
353	! be needed later on if reinitializing soil thermal inertia
354	surfithfi(i)=ithfi(i,1)
355	enddo
356	! also copy albedo and soil thermal inertia on (local) dynamics grid
357	! so that options below can manipulate either (but must then ensure
358	! to correctly recast things on physics grid)
359	call gr_fi_dyn(1,ngridmx,iip1,jjp1,albfi,alb)
360	call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
361	call gr_fi_dyn(1,ngridmx,iip1,jjp1,surfithfi,surfith)
362
363	endif
364	c-----------------------------------------------------------------------
365	c Initialisation des constantes dynamique
366	c-----------------------------------------------------------------------
367
368	! JVO 2017 : Added a +1 shift in tab_cntrl indexes to be coherent
369	! with Earth-like Titan start files ( cf dynetat0 from common )
370
371	kappa = tab_cntrl(10)
372	etot0 = tab_cntrl(13)
373	ptot0 = tab_cntrl(14)
374	ztot0 = tab_cntrl(15)
375	stot0 = tab_cntrl(16)
376	ang0 = tab_cntrl(17)
377	write(,) "Newstart: kappa,etot0,ptot0,ztot0,stot0,ang0"
378	write(,) kappa,etot0,ptot0,ztot0,stot0,ang0
379
380	! for vertical coordinate
381	preff=tab_cntrl(19) ! reference surface pressure
382	pa=tab_cntrl(18) ! reference pressure at which coord is purely pressure
383	!NB: in start_archive files tab_cntrl(17)=tab_cntrl(18)=0
384	write(,) "Newstart: preff=",preff," pa=",pa
385	yes=' '
386	do while ((yes.ne.'y').and.(yes.ne.'n'))
387	write(,) "Change the values of preff and pa? (y/n)"
388	read(*,fmt='(a)') yes
389	end do
390	if (yes.eq.'y') then
391	write(,)"New value of reference surface pressure preff?"
392	write(,)" (for Titan, typically preff=140000)"
393	read(,) preff
394	write(,)"New value of reference pressure pa for purely"
395	write(,)"pressure levels (for hybrid coordinates)?"
396	write(,)" (for Titan, typically pa=10000)"
397	read(,) pa
398	endif
399	c-----------------------------------------------------------------------
400	c Lecture du tab_cntrl et initialisation des constantes physiques
401	c - pour start: Lmodif = 0 => pas de modifications possibles
402	c (modif dans le tabfi de readfi + loin)
403	c - pour start_archive: Lmodif = 1 => modifications possibles
404	c-----------------------------------------------------------------------
405	if (choix_1.eq.0) then
406	! tabfi requires that input file be first opened by open_startphy(fichnom)
407	fichnom = 'start_archive.nc'
408	call open_startphy(fichnom)
409	call tabfi (ngridmx,nid,Lmodif,tab0,day_ini,lllm,p_rad,
410	. p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
411	else if (choix_1.eq.1) then
412	fichnom = 'startfi.nc'
413	call open_startphy(fichnom)
414	Lmodif=1 ! Lmodif set to 1 to allow modifications in phyeta0
415	call tabfi (ngridmx,nid_fi,Lmodif,tab0,day_ini,lllm,p_rad,
416	. p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
417	endif
418
419	if (p_omeg .eq. -9999.) then
420	p_omeg = 8.*atan(1.)/p_daysec
421	print*,"new value of omega",p_omeg
422	endif
423
424	rad = p_rad
425	omeg = p_omeg
426	g = p_g
427	cpp = p_cpp
428	mugaz = p_mugaz
429	daysec = p_daysec
430
431	if (p_omeg .eq. -9999.) then
432	p_omeg = 8.*atan(1.)/p_daysec
433	print*,"new value of omega",p_omeg
434	endif
435
436	kappa = 8.3141000./(p_mugazp_cpp) ! added by RDW
437
438	c=======================================================================
439	c INITIALISATIONS DIVERSES
440	c=======================================================================
441	! Load parameters from run.def file
442	CALL defrun_new( 99, .TRUE. )
443	! Initialize dynamics geometry and co. (which, when using start.nc may
444	! have changed because of modifications (tabi, preff, pa) above)
445	CALL iniconst
446	CALL inigeom
447	idum=-1
448	idum=0
449
450	! Initialize the physics for start_archive only
451	IF (choix_1.eq.0) THEN
452	CALL iniphysiq(iim,jjm,llm,
453	& (jjm-1)*iim+2,comm_lmdz,
454	& daysec,day_ini,dtphys,
455	& rlatu,rlatv,rlonu,rlonv,
456	& aire,cu,cv,rad,g,r,cpp,
457	& 1)
458	ENDIF
459
460	c=======================================================================
461	c lecture topographie, albedo, inertie thermique, relief sous-maille
462	c=======================================================================
463
464	if (choix_1.eq.0) then ! for start_archive files,
465	! where an external "surface.nc" file is needed
466
467	c do while((relief(1:3).ne.'mol').AND.(relief(1:3).ne.'pla'))
468	c write(,)
469	c write(,) 'choix du relief (mola,pla)'
470	c write(,) '(Topographie MGS MOLA, plat)'
471	c read(*,fmt='(a3)') relief
472	relief="mola"
473	c enddo
474
475	! First get the correct value of datadir, if not already done:
476	! default 'datadir' is set in "datafile_mod"
477	call getin_p("datadir",datadir)
478	write(,) 'Available surface data files are:'
479	filestring='ls '//trim(datadir)//'/'//
480	& trim(surfdir)//' \| grep .nc'
481	call system(filestring)
482	! but in ye old days these files were in datadir, so scan it as well
483	! for the sake of retro-compatibility
484	filestring='ls '//trim(datadir)//'/'//' \| grep .nc'
485	call system(filestring)
486
487	write(,) ''
488	write(,) 'Please choose the relevant file',
489	& ' (e.g. "surface_mars.nc")'
490	write(,) ' or "none" to not use any (and set planetary'
491	write(,) ' albedo and surface thermal inertia)'
492	read(*,fmt='(a50)') surfacefile
493
494	if (surfacefile.ne."none") then
495
496	CALL datareadnc(relief,surfacefile,phis,alb,surfith,
497	& zmeaS,zstdS,zsigS,zgamS,ztheS)
498	else
499	! specific case when not using a "surface.nc" file
500	phis(:,:)=0
501	zmeaS(:,:)=0
502	zstdS(:,:)=0
503	zsigS(:,:)=0
504	zgamS(:,:)=0
505	ztheS(:,:)=0
506
507	write(,) "Enter value of albedo of the bare ground:"
508	read(,) alb(1,1)
509	alb(:,:)=alb(1,1)
510
511	write(,) "Enter value of thermal inertia of soil:"
512	read(,) surfith(1,1)
513	surfith(:,:)=surfith(1,1)
514
515	endif ! of if (surfacefile.ne."none")
516
517	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
518	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,surfith,surfithfi)
519	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
520	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zmeaS,zmea)
521	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zstdS,zstd)
522	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zsigS,zsig)
523	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zgamS,zgam)
524	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ztheS,zthe)
525
526	endif ! of if (choix_1.eq.0)
527
528
529	c=======================================================================
530	c Lecture des fichiers (start ou start_archive)
531	c=======================================================================
532
533	if (choix_1.eq.0) then
534
535	write(,) 'Reading file START_ARCHIVE'
536	CALL lect_start_archive(ngridmx,llm,
537	& date,tsurf,tsoil,emis,q2,
538	& t,ucov,vcov,ps,teta,phisold_newgrid,
539	& q,qsurf,surfith,nid)
540	write(,) "OK, read start_archive file"
541	! copy soil thermal inertia
542	ithfi(:,:)=inertiedat(:,:)
543
544	ierr= NF_CLOSE(nid)
545
546	else if (choix_1.eq.1) then
547	!do nothing, start and startfi have already been read
548	else
549	CALL exit(1)
550	endif
551
552	dtvr = daysec/FLOAT(day_step)
553	dtphys = dtvr * FLOAT(iphysiq)
554
555	c=======================================================================
556	c
557	c=======================================================================
558
559	do ! infinite loop on list of changes
560
561	write(,)
562	write(,)
563	write(,) 'List of possible changes :'
564	write(,) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
565	write(,)
566	write(,) 'flat : no topography ("aquaplanet")'
567	write(,) 'set_ps_to_preff : used if changing preff with topo'
568	write(,) 'bilball : uniform albedo and thermal inertia'
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(,) 'isotherm : Isothermal Temperatures, wind set to zero'
575	write(,) 'radequi : Earth-like radiative equilibrium temperature
576	$ profile (lat-alt) and winds set to zero'
577	write(,) 'ptot : change total pressure'
578	write(,) 'emis : change surface emissivity'
579	write(,) 'therm_ini_s : Set soil thermal inertia to reference sur
580	&face values'
581	write(,) 'slab_ocean_0 : initialisation of slab
582	$ocean variables'
583
584	write(,)
585	write(,) 'Change to perform ?'
586	write(,) ' (enter keyword or return to end)'
587	write(,)
588
589	read(*,fmt='(a20)') modif
590	if (modif(1:1) .eq. ' ')then
591	exit ! exit loop on changes
592	endif
593
594	write(,)
595	write(,) trim(modif) , ' : '
596
597	c 'flat : no topography ("aquaplanet")'
598	c -------------------------------------
599	if (trim(modif) .eq. 'flat') then
600	c set topo to zero
601	z_reel(1:iip1,1:jjp1)=0
602	phis(1:iip1,1:jjp1)=g*z_reel(1:iip1,1:jjp1)
603	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
604	write(,) 'topography set to zero.'
605	write(,) 'WARNING : the subgrid topography parameters',
606	& ' were not set to zero ! => set calllott to F'
607
608	c Choice of surface pressure
609	yes=' '
610	do while ((yes.ne.'y').and.(yes.ne.'n'))
611	write(,) 'Do you wish to choose homogeneous surface',
612	& 'pressure (y) or let newstart interpolate ',
613	& ' the previous field (n)?'
614	read(*,fmt='(a)') yes
615	end do
616	if (yes.eq.'y') then
617	flagps0=.true.
618	write(,) 'New value for ps (Pa) ?'
619	201 read(,,iostat=ierr) patm
620	if(ierr.ne.0) goto 201
621	write(,) patm
622	if (patm.eq.-9999.) then
623	patm = preff
624	write(,) "we set patm = preff", preff
625	endif
626	write(,)
627	write(,) ' new ps everywhere (Pa) = ', patm
628	write(,)
629	do j=1,jjp1
630	do i=1,iip1
631	ps(i,j)=patm
632	enddo
633	enddo
634	end if
635
636	c 'set_ps_to_preff' : used if changing preff with topo
637	c ----------------------------------------------------
638	else if (trim(modif) .eq. 'set_ps_to_preff') then
639	do j=1,jjp1
640	do i=1,iip1
641	ps(i,j)=preff
642	enddo
643	enddo
644
645
646	c bilball : uniform albedo, thermal inertia
647	c -----------------------------------------
648	else if (trim(modif) .eq. 'bilball') then
649	write(,) 'constante albedo and iner.therm:'
650	write(,) 'New value for albedo (ex: 0.25) ?'
651	101 read(,,iostat=ierr) alb_bb
652	if(ierr.ne.0) goto 101
653	write(,)
654	write(,) ' uniform albedo (new value):',alb_bb
655	write(,)
656
657	write(,) 'New value for thermal inertia (eg: 247) ?'
658	102 read(,,iostat=ierr) ith_bb
659	if(ierr.ne.0) goto 102
660	write(,) 'uniform thermal inertia (new value):',ith_bb
661	DO j=1,jjp1
662	DO i=1,iip1
663	alb(i,j) = alb_bb ! albedo
664	do isoil=1,nsoilmx
665	ith(i,j,isoil) = ith_bb ! thermal inertia
666	enddo
667	END DO
668	END DO
669	! CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ith,ithfi)
670	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
671	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
672
673	c ptot : Modification of the total pressure: ice + current atmosphere
674	c -------------------------------------------------------------------
675	else if (trim(modif).eq.'ptot') then
676
677	c calcul de la pression totale glace + atm actuelle
678	patm=0.
679	airetot=0.
680	pcap=0.
681	DO j=1,jjp1
682	DO i=1,iim
683	ig=1+(j-2)*iim +i
684	if(j.eq.1) ig=1
685	if(j.eq.jjp1) ig=ngridmx
686	patm = patm + ps(i,j)*aire(i,j)
687	airetot= airetot + aire(i,j)
688	ENDDO
689	ENDDO
690	ptoto = pcap + patm
691
692	print*,'Current total pressure at surface (atm) ',
693	& ptoto/airetot
694
695	print*,'new value?'
696	read(,) ptotn
697	ptotn=ptotn*airetot
698	patmn=ptotn-pcap
699	print*,'ptoto,patm,ptotn,patmn'
700	print*,ptoto,patm,ptotn,patmn
701	print*,'Mult. factor for pressure (atm only)', patmn/patm
702	do j=1,jjp1
703	do i=1,iip1
704	ps(i,j)=ps(i,j)*patmn/patm
705	enddo
706	enddo
707
708
709
710	c Correction pour la conservation des traceurs
711	yes=' '
712	do while ((yes.ne.'y').and.(yes.ne.'n'))
713	write(,) 'Do you wish to conserve tracer total mass (y)',
714	& ' or tracer mixing ratio (n) ?'
715	read(*,fmt='(a)') yes
716	end do
717
718	if (yes.eq.'y') then
719	write(,) 'OK : conservation of tracer total mass'
720	DO iq =1, nqtot
721	DO l=1,llm
722	DO j=1,jjp1
723	DO i=1,iip1
724	q(i,j,l,iq)=q(i,j,l,iq)*patm/patmn
725	ENDDO
726	ENDDO
727	ENDDO
728	ENDDO
729	else
730	write(,) 'OK : conservation of tracer mixing ratio'
731	end if
732
733	c Correction pour la pression au niveau de la mer
734	yes=' '
735	do while ((yes.ne.'y').and.(yes.ne.'n'))
736	write(,) 'Do you wish fix pressure at sea level (y)',
737	& ' or not (n) ?'
738	read(*,fmt='(a)') yes
739	end do
740
741	if (yes.eq.'y') then
742	write(,) 'Value?'
743	read(,,iostat=ierr) psea
744	DO i=1,iip1
745	DO j=1,jjp1
746	ps(i,j)=psea
747
748	ENDDO
749	ENDDO
750	write(,) 'psea=',psea
751	else
752	write(,) 'no'
753	end if
754
755
756	c emis : change surface emissivity (added by RW)
757	c ----------------------------------------------
758	else if (trim(modif).eq.'emis') then
759
760	print*,'new value?'
761	read(,) emisread
762
763	do i=1,ngridmx
764	emis(i)=emisread
765	enddo
766
767	c qname : change tracer name
768	c --------------------------
769	else if (trim(modif).eq.'qname') then
770	yes='y'
771	do while (yes.eq.'y')
772	write(,) 'Which tracer name do you want to change ?'
773	do iq=1,nqtot
774	write(*,'(i3,a3,a20)')iq,' : ',trim(tname(iq))
775	enddo
776	write(*,'(a35,i3)')
777	& '(enter tracer number; between 1 and ',nqtot
778	write(,)' or any other value to quit this option)'
779	read(,) iq
780	if ((iq.ge.1).and.(iq.le.nqtot)) then
781	write(,)'Change tracer name ',trim(tname(iq)),' to ?'
782	read(,) txt
783	tname(iq)=txt
784	write(,)'Do you want to change another tracer name (y/n)?'
785	read(*,'(a)') yes
786	else
787	! inapropiate value of iq; quit this option
788	yes='n'
789	endif ! of if ((iq.ge.1).and.(iq.le.nqtot))
790	enddo ! of do while (yes.ne.'y')
791
792	c q=0 : set tracers to zero
793	c -------------------------
794	else if (trim(modif).eq.'q=0') then
795	c mise a 0 des q (traceurs)
796	write(,) 'Tracers set to 0 (1.E-30 in fact)'
797	DO iq =1, nqtot
798	DO l=1,llm
799	DO j=1,jjp1
800	DO i=1,iip1
801	q(i,j,l,iq)=1.e-30
802	ENDDO
803	ENDDO
804	ENDDO
805	ENDDO
806
807	c set surface tracers to zero
808	DO iq =1, nqtot
809	DO ig=1,ngridmx
810	qsurf(ig,iq)=0.
811	ENDDO
812	ENDDO
813
814	c q=x : initialise tracer manually
815	c --------------------------------
816	else if (trim(modif).eq.'q=x') then
817	write(,) 'Which tracer do you want to modify ?'
818	do iq=1,nqtot
819	write(,)iq,' : ',trim(tname(iq))
820	enddo
821	write(,) '(choose between 1 and ',nqtot,')'
822	read(,) iq
823	write(,)'mixing ratio of tracer ',trim(tname(iq)),
824	& ' ? (kg/kg)'
825	read(,) val
826	DO l=1,llm
827	DO j=1,jjp1
828	DO i=1,iip1
829	q(i,j,l,iq)=val
830	ENDDO
831	ENDDO
832	ENDDO
833	write(,) 'SURFACE value of tracer ',trim(tname(iq)),
834	& ' ? (kg/m2)'
835	read(,) val
836	DO ig=1,ngridmx
837	qsurf(ig,iq)=val
838	ENDDO
839
840	c t=profile : initialize temperature with a given profile
841	c -------------------------------------------------------
842	else if (trim(modif) .eq. 't=profile') then
843	write(,) 'Temperature profile from ASCII file'
844	write(,) "'profile.in' e.g. 1D output"
845	write(,) "(one value per line in file; starting with"
846	write(,) "surface value, the 1st atmospheric layer"
847	write(,) "followed by 2nd, etc. up to top of atmosphere)"
848	txt="profile.in"
849	open(33,file=trim(txt),status='old',form='formatted',
850	& iostat=ierr)
851	if (ierr.eq.0) then
852	! OK, found file 'profile_...', load the profile
853	do l=1,llm+1
854	read(33,*,iostat=ierr) profile(l)
855	write(,) profile(l)
856	if (ierr.ne.0) then ! something went wrong
857	exit ! quit loop
858	endif
859	enddo
860	if (ierr.eq.0) then
861	tsurf(1:ngridmx)=profile(1)
862	tsoil(1:ngridmx,1:nsoilmx)=profile(1)
863	do l=1,llm
864	Tset(1:iip1,1:jjp1,l)=profile(l+1)
865	flagtset=.true.
866	enddo
867	ucov(1:iip1,1:jjp1,1:llm)=0.
868	vcov(1:iip1,1:jjm,1:llm)=0.
869	q2(1:ngridmx,1:llm+1)=0.
870	else
871	write(,)'problem reading file ',trim(txt),' !'
872	write(,)'No modifications to temperature'
873	endif
874	else
875	write(,)'Could not find file ',trim(txt),' !'
876	write(,)'No modifications to temperature'
877	endif
878
879	c q=profile : initialize tracer with a given profile
880	c --------------------------------------------------
881	else if (trim(modif) .eq. 'q=profile') then
882	write(,) 'Tracer profile will be sought in ASCII file'
883	write(,) "'profile_tracer' where 'tracer' is tracer name"
884	write(,) "(one value per line in file; starting with"
885	write(,) "surface value, the 1st atmospheric layer"
886	write(,) "followed by 2nd, etc. up to top of atmosphere)"
887	write(,) 'Which tracer do you want to set?'
888	do iq=1,nqtot
889	write(,)iq,' : ',trim(tname(iq))
890	enddo
891	write(,) '(choose between 1 and ',nqtot,')'
892	read(,) iq
893	if ((iq.lt.1).or.(iq.gt.nqtot)) then
894	! wrong value for iq, go back to menu
895	write(,) "wrong input value:",iq
896	cycle
897	endif
898	! look for input file 'profile_tracer'
899	txt="profile_"//trim(tname(iq))
900	open(41,file=trim(txt),status='old',form='formatted',
901	& iostat=ierr)
902	if (ierr.eq.0) then
903	! OK, found file 'profile_...', load the profile
904	do l=1,llm+1
905	read(41,*,iostat=ierr) profile(l)
906	if (ierr.ne.0) then ! something went wrong
907	exit ! quit loop
908	endif
909	enddo
910	if (ierr.eq.0) then
911	! initialize tracer values
912	qsurf(:,iq)=profile(1)
913	do l=1,llm
914	q(:,:,l,iq)=profile(l+1)
915	enddo
916	write(,)'OK, tracer ',trim(tname(iq)),' initialized '
917	& ,'using values from file ',trim(txt)
918	else
919	write(,)'problem reading file ',trim(txt),' !'
920	write(,)'No modifications to tracer ',trim(tname(iq))
921	endif
922	else
923	write(,)'Could not find file ',trim(txt),' !'
924	write(,)'No modifications to tracer ',trim(tname(iq))
925	endif
926
927
928	c isotherm : Isothermal temperatures and no winds
929	c -----------------------------------------------
930	else if (trim(modif) .eq. 'isotherm') then
931
932	write(,)'Isothermal temperature of the atmosphere,
933	& surface and subsurface'
934	write(,) 'Value of this temperature ? :'
935	203 read(,,iostat=ierr) Tiso
936	if(ierr.ne.0) goto 203
937
938	tsurf(1:ngridmx)=Tiso
939
940	tsoil(1:ngridmx,1:nsoilmx)=Tiso
941
942	Tset(1:iip1,1:jjp1,1:llm)=Tiso
943	flagtset=.true.
944
945	t(1:iip1,1:jjp1,1:llm)=Tiso
946	!! otherwise hydrost. integrations below
947	!! use the wrong temperature
948	!! -- NB: Tset might be useless
949
950	ucov(1:iip1,1:jjp1,1:llm)=0
951	vcov(1:iip1,1:jjm,1:llm)=0
952	q2(1:ngridmx,1:llm+1)=0
953
954	c radequi : Radiative equilibrium profile of temperatures and no winds
955	c --------------------------------------------------------------------
956	else if (trim(modif) .eq. 'radequi') then
957
958	write(,)'radiative equilibrium temperature profile'
959
960	do ig=1, ngridmx
961	teque= 335.0-60.0sin(latfi(ig))sin(latfi(ig))-
962	& 10.0cos(latfi(ig))cos(latfi(ig))
963	tsurf(ig) = MAX(220.0,teque)
964	end do
965	do l=2,nsoilmx
966	do ig=1, ngridmx
967	tsoil(ig,l) = tsurf(ig)
968	end do
969	end do
970	DO j=1,jjp1
971	DO i=1,iim
972	Do l=1,llm
973	teque=335.-60.0sin(rlatu(j))sin(rlatu(j))-
974	& 10.0cos(rlatu(j))cos(rlatu(j))
975	Tset(i,j,l)=MAX(220.0,teque)
976	end do
977	end do
978	end do
979	flagtset=.true.
980	ucov(1:iip1,1:jjp1,1:llm)=0
981	vcov(1:iip1,1:jjm,1:llm)=0
982	q2(1:ngridmx,1:llm+1)=0
983
984	! therm_ini_s: (re)-set soil thermal inertia to reference surface values
985	! ----------------------------------------------------------------------
986
987	else if (trim(modif) .eq. 'therm_ini_s') then
988	! write(,)"surfithfi(1):",surfithfi(1)
989	do isoil=1,nsoilmx
990	inertiedat(1:ngridmx,isoil)=surfithfi(1:ngridmx)
991	enddo
992	write(,)'OK: Soil thermal inertia has been reset to referenc
993	&e surface values'
994	! write(,)"inertiedat(1,1):",inertiedat(1,1)
995	ithfi(:,:)=inertiedat(:,:)
996	! recast ithfi() onto ith()
997	call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
998	! Check:
999	! do i=1,iip1
1000	! do j=1,jjp1
1001	! do isoil=1,nsoilmx
1002	! write(77,*) i,j,isoil," ",ith(i,j,isoil)
1003	! enddo
1004	! enddo
1005	! enddo
1006
1007
1008
1009	c slab_ocean_initialisation
1010	c ------------------------------------------------
1011	else if (trim(modif) .eq. 'slab_ocean_0') then
1012
1013	write(,) "No ocean yet on Titan ! Can't use this option"
1014	stop
1015
1016	else
1017	write(,) ' Unknown (misspelled?) option!!!'
1018	end if ! of if (trim(modif) .eq. '...') elseif ...
1019
1020
1021
1022	enddo ! of do ! infinite loop on liste of changes
1023
1024	999 continue
1025
1026
1027	c================================================================
1028	c Initialisation for methane surface tank (added by JVO 2017)
1029	c================================================================
1030	DO ig=1, ngridmx
1031	tankCH4(ig)=2.0
1032	ENDDO
1033
1034	c=======================================================================
1035	c Correct pressure on the new grid (menu 0)
1036	c=======================================================================
1037
1038
1039	if ((choix_1.eq.0).and.(.not.flagps0)) then
1040	r = 1000.*8.31/mugaz
1041
1042	do j=1,jjp1
1043	do i=1,iip1
1044	ps(i,j) = ps(i,j) *
1045	. exp((phisold_newgrid(i,j)-phis(i,j)) /
1046	. (t(i,j,1) * r))
1047	end do
1048	end do
1049
1050	c periodicite de ps en longitude
1051	do j=1,jjp1
1052	ps(1,j) = ps(iip1,j)
1053	end do
1054	end if
1055
1056
1057	c=======================================================================
1058	c=======================================================================
1059
1060	c=======================================================================
1061	c Initialisation de la physique / ecriture de newstartfi :
1062	c=======================================================================
1063
1064
1065	CALL inifilr
1066	CALL pression(ip1jmp1, ap, bp, ps, p3d)
1067
1068
1069	c=========================================================================
1070	c Calcul de la dimension verticale pour la chimie - JVO 2017
1071	c start_archive seulement, la grille verticale pouvant avoir ete modifiee
1072	c==========================================================================
1073
1074	IF (choix_1.eq.0) THEN
1075
1076	! Calculate the # of upper chemistry layers with the "new" pressure grid
1077	! For this we use Vervack profile for upper atmosphere with dz=10km
1078
1079	CALL gr_kim_vervack
1080
1081	WRITE(,)
1082	WRITE(,) " With the compiled vertical grid we found :"
1083	WRITE(,) " Number of upper chemistry layers =", nlaykim_up
1084
1085	! Regriding is then done, if needed
1086
1087	IF (nlaykimold.ne.nlaykim_up) THEN
1088
1089	WRITE(,) " Warning, nlaykimold=", nlaykimold
1090	WRITE(,) ' which implies that a regriding on upper chemistry
1091	& will be performed.'
1092	WRITE(,)
1093
1094	! CALL regrid_kim
1095
1096	ENDIF
1097
1098	endif ! of if (choix_1.eq.0)
1099
1100
1101	c-----------------------------------------------------------------------
1102	c Initialisation pks:
1103	c-----------------------------------------------------------------------
1104
1105	CALL exner_hyb(ip1jmp1, ps, p3d, beta, pks, pk, pkf)
1106	! Calcul de la temperature potentielle teta
1107
1108	if (flagtset) then
1109	DO l=1,llm
1110	DO j=1,jjp1
1111	DO i=1,iim
1112	teta(i,j,l) = Tset(i,j,l) * cpp/pk(i,j,l)
1113	ENDDO
1114	teta (iip1,j,l)= teta (1,j,l)
1115	ENDDO
1116	ENDDO
1117	else if (choix_1.eq.0) then
1118	DO l=1,llm
1119	DO j=1,jjp1
1120	DO i=1,iim
1121	teta(i,j,l) = t(i,j,l) * cpp/pk(i,j,l)
1122	ENDDO
1123	teta (iip1,j,l)= teta (1,j,l)
1124	ENDDO
1125	ENDDO
1126	endif
1127
1128	C Calcul intermediaire
1129	c
1130	if (choix_1.eq.0) then
1131	CALL massdair( p3d, masse )
1132	c
1133	print *,' ALPHAX ',alphax
1134
1135	DO l = 1, llm
1136	DO i = 1, iim
1137	xppn(i) = aire( i, 1 ) * masse( i , 1 , l )
1138	xpps(i) = aire( i,jjp1 ) * masse( i , jjp1 , l )
1139	ENDDO
1140	xpn = SUM(xppn)/apoln
1141	xps = SUM(xpps)/apols
1142	DO i = 1, iip1
1143	masse( i , 1 , l ) = xpn
1144	masse( i , jjp1 , l ) = xps
1145	ENDDO
1146	ENDDO
1147	endif
1148	phis(iip1,:) = phis(1,:)
1149
1150	itau=0
1151	if (choix_1.eq.0) then
1152	day_ini=int(date)
1153	endif
1154
1155	c
1156	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
1157
1158	CALL caldyn0( itau,ucov,vcov,teta,ps,masse,pk,phis ,
1159	* phi,w, pbaru,pbarv,day_ini+time )
1160
1161
1162	CALL dynredem0("restart.nc",day_ini,phis)
1163	CALL dynredem1("restart.nc",0.0,vcov,ucov,teta,q,masse,ps)
1164	C
1165	C Ecriture etat initial physique
1166	C
1167
1168	call physdem0("restartfi.nc",lonfi,latfi,nsoilmx,ngridmx,llm,
1169	& nqtot,dtphys,real(day_ini),0.0,
1170	& airefi,albfi,ithfi,zmea,zstd,zsig,zgam,zthe)
1171	call physdem1("restartfi.nc",nsoilmx,ngridmx,llm,nqtot,
1172	& dtphys,real(day_ini),
1173	& tsurf,tsoil,emis,q2,qsurf,tankCH4)
1174
1175	c=======================================================================
1176	c Formats
1177	c=======================================================================
1178
1179	1 FORMAT(//10x,'la valeur de im =',i4,2x,'lue sur le fichier de dema
1180	*rrage est differente de la valeur parametree iim =',i4//)
1181	2 FORMAT(//10x,'la valeur de jm =',i4,2x,'lue sur le fichier de dema
1182	*rrage est differente de la valeur parametree jjm =',i4//)
1183	3 FORMAT(//10x,'la valeur de lllm =',i4,2x,'lue sur le fichier demar
1184	*rage est differente de la valeur parametree llm =',i4//)
1185
1186	write(,) "newstart: All is well that ends well."
1187
1188	end
1189

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

Download in other formats:

Original Format