Context Navigation

newstart.F @ 1891

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

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

Download in other formats:

Original Format