Context Navigation

newstart.F @ 3595

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

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

Download in other formats:

Original Format