Context Navigation

newstart.F @ 3026

Last change on this file since 3026 was 2366, checked in by jvatant, 5 years ago

Titan GCM : Major maintenance catching up commits from the generic including :

r2356 and 2354 removing obsolete old dynamical core
various minor addition to physics and gestion of phys_state_var_mode, especially in dyn1d
adding MESOSCALE CPP keys around chemistry and microphysics (disabled in mesoscale for now)

File size: 41.1 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
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	implicit none
44
45	include "dimensions.h"
46	integer, parameter :: ngridmx = (2+(jjm-1)*iim - 1/jjm)
47	include "paramet.h"
48	include "comgeom2.h"
49	include "comdissnew.h"
50	include "netcdf.inc"
51
52	c=======================================================================
53	c Declarations
54	c=======================================================================
55
56	c Variables dimension du fichier "start_archive"
57	c------------------------------------
58	CHARACTER relief*3
59
60
61	c Variables pour les lectures NetCDF des fichiers "start_archive"
62	c--------------------------------------------------
63	INTEGER nid_dyn, nid_fi,nid,nvarid
64	INTEGER length
65	parameter (length = 100)
66	INTEGER tab0
67	INTEGER NB_ETATMAX
68	parameter (NB_ETATMAX = 100)
69
70	REAL date
71	REAL p_rad,p_omeg,p_g,p_cpp,p_mugaz,p_daysec
72
73	c Variable histoire
74	c------------------
75	REAL vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm) ! vents covariants
76	REAL phis(iip1,jjp1)
77	REAL,ALLOCATABLE :: q(:,:,:,:) ! champs advectes
78
79	c autre variables dynamique nouvelle grille
80	c------------------------------------------
81	REAL pks(iip1,jjp1)
82	REAL w(iip1,jjp1,llm+1)
83	REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
84	! REAL dv(ip1jm,llm),du(ip1jmp1,llm)
85	! REAL dh(ip1jmp1,llm),dp(ip1jmp1)
86	REAL phi(iip1,jjp1,llm)
87
88	integer klatdat,klongdat
89	PARAMETER (klatdat=180,klongdat=360)
90
91	c Physique sur grille scalaire
92	c----------------------------
93	real zmeaS(iip1,jjp1),zstdS(iip1,jjp1)
94	real zsigS(iip1,jjp1),zgamS(iip1,jjp1),ztheS(iip1,jjp1)
95
96	c variable physique
97	c------------------
98	REAL tsurf(ngridmx) ! surface temperature
99	REAL,ALLOCATABLE :: tsoil(:,:) ! soil temperature
100	REAL emis(ngridmx) ! surface emissivity
101	real emisread ! added by RW
102	REAL,ALLOCATABLE :: qsurf(:,:)
103	REAL q2(ngridmx,llm+1)
104	! REAL rnaturfi(ngridmx)
105	real alb(iip1,jjp1),albfi(ngridmx) ! albedos
106	real,ALLOCATABLE :: ith(:,:,:),ithfi(:,:) ! thermal inertia (3D)
107	real surfith(iip1,jjp1),surfithfi(ngridmx) ! surface thermal inertia (2D)
108	REAL latfi(ngridmx),lonfi(ngridmx),airefi(ngridmx)
109
110	INTEGER i,j,l,isoil,ig,idum
111	real mugaz ! molar mass of the atmosphere
112
113	integer ierr
114
115	c Variables on the new grid along scalar points
116	c------------------------------------------------------
117	! REAL p(iip1,jjp1)
118	REAL t(iip1,jjp1,llm)
119	REAL tset(iip1,jjp1,llm)
120	real phisold_newgrid(iip1,jjp1)
121	REAL :: teta(iip1, jjp1, llm)
122	REAL :: pk(iip1,jjp1,llm)
123	REAL :: pkf(iip1,jjp1,llm)
124	REAL :: ps(iip1, jjp1)
125	REAL :: masse(iip1,jjp1,llm)
126	REAL :: xpn,xps,xppn(iim),xpps(iim)
127	REAL :: p3d(iip1, jjp1, llm+1)
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	! Take care of arrays in common modules
269	! ALLOCATE ARRAYS in surfdat_h (if not already done, e.g. when using start_archive)
270	IF (.not. ALLOCATED(albedodat)) ALLOCATE(albedodat(ngridmx))
271	IF (.not. ALLOCATED(phisfi)) ALLOCATE(phisfi(ngridmx))
272	IF (.not. ALLOCATED(zmea)) ALLOCATE(zmea(ngridmx))
273	IF (.not. ALLOCATED(zstd)) ALLOCATE(zstd(ngridmx))
274	IF (.not. ALLOCATED(zsig)) ALLOCATE(zsig(ngridmx))
275	IF (.not. ALLOCATED(zgam)) ALLOCATE(zgam(ngridmx))
276	IF (.not. ALLOCATED(zthe)) ALLOCATE(zthe(ngridmx))
277
278	c-----------------------------------------------------------------------
279	c Lecture du tableau des parametres du run (pour la dynamique)
280	c-----------------------------------------------------------------------
281	if (choix_1.eq.0) then
282
283	write(,) 'reading tab_cntrl START_ARCHIVE'
284	c
285	ierr = NF_INQ_VARID (nid, "controle", nvarid)
286	#ifdef NC_DOUBLE
287	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tab_cntrl)
288	#else
289	ierr = NF_GET_VAR_REAL(nid, nvarid, tab_cntrl)
290	#endif
291	c
292	else if (choix_1.eq.1) then
293
294	write(,) 'reading tab_cntrl START'
295	c
296	ierr = NF_INQ_VARID (nid_dyn, "controle", nvarid)
297	#ifdef NC_DOUBLE
298	ierr = NF_GET_VAR_DOUBLE(nid_dyn, nvarid, tab_cntrl)
299	#else
300	ierr = NF_GET_VAR_REAL(nid_dyn, nvarid, tab_cntrl)
301	#endif
302	c
303	write(,) 'reading tab_cntrl STARTFI'
304	c
305	ierr = NF_INQ_VARID (nid_fi, "controle", nvarid)
306	#ifdef NC_DOUBLE
307	ierr = NF_GET_VAR_DOUBLE(nid_fi, nvarid, tab_cntrl_bis)
308	#else
309	ierr = NF_GET_VAR_REAL(nid_fi, nvarid, tab_cntrl_bis)
310	#endif
311	c
312	do i=1,50
313	tab_cntrl(i+50)=tab_cntrl_bis(i)
314	enddo
315	write(,) 'printing tab_cntrl', tab_cntrl
316	do i=1,100
317	write(,) i,tab_cntrl(i)
318	enddo
319
320	write(,) 'Reading file START'
321	fichnom = 'start.nc'
322	CALL dynetat0(fichnom,vcov,ucov,teta,q,masse,
323	. ps,phis,time)
324
325	CALL iniconst
326	CALL inigeom
327
328	! Initialize the physics
329	CALL iniphysiq(iim,jjm,llm,
330	& (jjm-1)*iim+2,comm_lmdz,
331	& daysec,day_ini,dtphys,
332	& rlatu,rlatv,rlonu,rlonv,
333	& aire,cu,cv,rad,g,r,cpp,
334	& 1)
335
336	! Initialize global tracer indexes (stored in tracer.h)
337	! ... this has to be done before phyetat0
338	CALL initracer2(nqtot,tname)
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	! Initialize global tracer indexes (stored in tracer.h)
461	CALL initracer2(nqtot,tname)
462
463	ENDIF
464
465	c=======================================================================
466	c lecture topographie, albedo, inertie thermique, relief sous-maille
467	c=======================================================================
468
469	if (choix_1.eq.0) then ! for start_archive files,
470	! where an external "surface.nc" file is needed
471
472	c do while((relief(1:3).ne.'mol').AND.(relief(1:3).ne.'pla'))
473	c write(,)
474	c write(,) 'choix du relief (mola,pla)'
475	c write(,) '(Topographie MGS MOLA, plat)'
476	c read(*,fmt='(a3)') relief
477	relief="mola"
478	c enddo
479
480	! First get the correct value of datadir, if not already done:
481	! default 'datadir' is set in "datafile_mod"
482	call getin_p("datadir",datadir)
483	write(,) 'Available surface data files are:'
484	filestring='ls '//trim(datadir)//'/'//
485	& trim(surfdir)//' \| grep .nc'
486	call system(filestring)
487	! but in ye old days these files were in datadir, so scan it as well
488	! for the sake of retro-compatibility
489	filestring='ls '//trim(datadir)//'/'//' \| grep .nc'
490	call system(filestring)
491
492	write(,) ''
493	write(,) 'Please choose the relevant file',
494	& ' (e.g. "surface_mars.nc")'
495	write(,) ' or "none" to not use any (and set planetary'
496	write(,) ' albedo and surface thermal inertia)'
497	read(*,fmt='(a50)') surfacefile
498
499	if (surfacefile.ne."none") then
500
501	CALL datareadnc(relief,surfacefile,phis,alb,surfith,
502	& zmeaS,zstdS,zsigS,zgamS,ztheS)
503	else
504	! specific case when not using a "surface.nc" file
505	phis(:,:)=0
506	zmeaS(:,:)=0
507	zstdS(:,:)=0
508	zsigS(:,:)=0
509	zgamS(:,:)=0
510	ztheS(:,:)=0
511
512	write(,) "Enter value of albedo of the bare ground:"
513	read(,) alb(1,1)
514	alb(:,:)=alb(1,1)
515
516	write(,) "Enter value of thermal inertia of soil:"
517	read(,) surfith(1,1)
518	surfith(:,:)=surfith(1,1)
519
520	endif ! of if (surfacefile.ne."none")
521
522	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
523	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,surfith,surfithfi)
524	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
525	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zmeaS,zmea)
526	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zstdS,zstd)
527	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zsigS,zsig)
528	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zgamS,zgam)
529	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ztheS,zthe)
530
531	endif ! of if (choix_1.eq.0)
532
533
534	c=======================================================================
535	c Lecture des fichiers (start ou start_archive)
536	c=======================================================================
537
538	if (choix_1.eq.0) then
539
540	write(,) 'Reading file START_ARCHIVE'
541	CALL lect_start_archive(ngridmx,llm,
542	& date,tsurf,tsoil,emis,q2,
543	& t,ucov,vcov,ps,teta,phisold_newgrid,
544	& q,qsurf,tankCH4,surfith,nid)
545	write(,) "OK, read start_archive file"
546	! copy soil thermal inertia
547	ithfi(:,:)=inertiedat(:,:)
548
549	ierr= NF_CLOSE(nid)
550
551	else if (choix_1.eq.1) then
552	!do nothing, start and startfi have already been read
553	else
554	CALL exit(1)
555	endif
556
557	dtvr = daysec/FLOAT(day_step)
558	dtphys = dtvr * FLOAT(iphysiq)
559
560	c=======================================================================
561	c
562	c=======================================================================
563
564	do ! infinite loop on list of changes
565
566	write(,)
567	write(,)
568	write(,) 'List of possible changes :'
569	write(,) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
570	write(,)
571	write(,) 'flat : no topography ("aquaplanet")'
572	write(,) 'set_ps_to_preff : used if changing preff with topo'
573	write(,) 'bilball : uniform albedo and thermal inertia'
574	write(,) 'qname : change tracer name'
575	write(,) 't=profile : read temperature profile in profile.in'
576	write(,) 'q=0 : ALL tracer =zero'
577	write(,) 'q=x : give a specific uniform value to one tracer'
578	write(,) 'q=profile : specify a profile for a tracer'
579	write(,) 'isotherm : Isothermal Temperatures, wind set to zero'
580	write(,) 'radequi : Earth-like radiative equilibrium temperature
581	$ profile (lat-alt) and winds set to zero'
582	write(,) 'ptot : change total pressure'
583	write(,) 'emis : change surface emissivity'
584	write(,) 'therm_ini_s : Set soil thermal inertia to reference sur
585	&face values'
586	write(,) 'slab_ocean_0 : initialisation of slab
587	$ocean variables'
588
589	write(,)
590	write(,) 'Change to perform ?'
591	write(,) ' (enter keyword or return to end)'
592	write(,)
593
594	read(*,fmt='(a20)') modif
595	if (modif(1:1) .eq. ' ')then
596	exit ! exit loop on changes
597	endif
598
599	write(,)
600	write(,) trim(modif) , ' : '
601
602	c 'flat : no topography ("aquaplanet")'
603	c -------------------------------------
604	if (trim(modif) .eq. 'flat') then
605	c set topo to zero
606	z_reel(1:iip1,1:jjp1)=0
607	phis(1:iip1,1:jjp1)=g*z_reel(1:iip1,1:jjp1)
608	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
609	write(,) 'topography set to zero.'
610	write(,) 'WARNING : the subgrid topography parameters',
611	& ' were not set to zero ! => set calllott to F'
612
613	c Choice of surface pressure
614	yes=' '
615	do while ((yes.ne.'y').and.(yes.ne.'n'))
616	write(,) 'Do you wish to choose homogeneous surface',
617	& 'pressure (y) or let newstart interpolate ',
618	& ' the previous field (n)?'
619	read(*,fmt='(a)') yes
620	end do
621	if (yes.eq.'y') then
622	flagps0=.true.
623	write(,) 'New value for ps (Pa) ?'
624	201 read(,,iostat=ierr) patm
625	if(ierr.ne.0) goto 201
626	write(,) patm
627	if (patm.eq.-9999.) then
628	patm = preff
629	write(,) "we set patm = preff", preff
630	endif
631	write(,)
632	write(,) ' new ps everywhere (Pa) = ', patm
633	write(,)
634	do j=1,jjp1
635	do i=1,iip1
636	ps(i,j)=patm
637	enddo
638	enddo
639	end if
640
641	c 'set_ps_to_preff' : used if changing preff with topo
642	c ----------------------------------------------------
643	else if (trim(modif) .eq. 'set_ps_to_preff') then
644	do j=1,jjp1
645	do i=1,iip1
646	ps(i,j)=preff
647	enddo
648	enddo
649
650
651	c bilball : uniform albedo, thermal inertia
652	c -----------------------------------------
653	else if (trim(modif) .eq. 'bilball') then
654	write(,) 'constante albedo and iner.therm:'
655	write(,) 'New value for albedo (ex: 0.25) ?'
656	101 read(,,iostat=ierr) alb_bb
657	if(ierr.ne.0) goto 101
658	write(,)
659	write(,) ' uniform albedo (new value):',alb_bb
660	write(,)
661
662	write(,) 'New value for thermal inertia (eg: 247) ?'
663	102 read(,,iostat=ierr) ith_bb
664	if(ierr.ne.0) goto 102
665	write(,) 'uniform thermal inertia (new value):',ith_bb
666	DO j=1,jjp1
667	DO i=1,iip1
668	alb(i,j) = alb_bb ! albedo
669	do isoil=1,nsoilmx
670	ith(i,j,isoil) = ith_bb ! thermal inertia
671	enddo
672	END DO
673	END DO
674	! CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ith,ithfi)
675	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
676	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
677
678	c ptot : Modification of the total pressure: ice + current atmosphere
679	c -------------------------------------------------------------------
680	else if (trim(modif).eq.'ptot') then
681
682	c calcul de la pression totale glace + atm actuelle
683	patm=0.
684	airetot=0.
685	pcap=0.
686	DO j=1,jjp1
687	DO i=1,iim
688	ig=1+(j-2)*iim +i
689	if(j.eq.1) ig=1
690	if(j.eq.jjp1) ig=ngridmx
691	patm = patm + ps(i,j)*aire(i,j)
692	airetot= airetot + aire(i,j)
693	ENDDO
694	ENDDO
695	ptoto = pcap + patm
696
697	print*,'Current total pressure at surface (atm) ',
698	& ptoto/airetot
699
700	print*,'new value?'
701	read(,) ptotn
702	ptotn=ptotn*airetot
703	patmn=ptotn-pcap
704	print*,'ptoto,patm,ptotn,patmn'
705	print*,ptoto,patm,ptotn,patmn
706	print*,'Mult. factor for pressure (atm only)', patmn/patm
707	do j=1,jjp1
708	do i=1,iip1
709	ps(i,j)=ps(i,j)*patmn/patm
710	enddo
711	enddo
712
713
714
715	c Correction pour la conservation des traceurs
716	yes=' '
717	do while ((yes.ne.'y').and.(yes.ne.'n'))
718	write(,) 'Do you wish to conserve tracer total mass (y)',
719	& ' or tracer mixing ratio (n) ?'
720	read(*,fmt='(a)') yes
721	end do
722
723	if (yes.eq.'y') then
724	write(,) 'OK : conservation of tracer total mass'
725	DO iq =1, nqtot
726	DO l=1,llm
727	DO j=1,jjp1
728	DO i=1,iip1
729	q(i,j,l,iq)=q(i,j,l,iq)*patm/patmn
730	ENDDO
731	ENDDO
732	ENDDO
733	ENDDO
734	else
735	write(,) 'OK : conservation of tracer mixing ratio'
736	end if
737
738	c Correction pour la pression au niveau de la mer
739	yes=' '
740	do while ((yes.ne.'y').and.(yes.ne.'n'))
741	write(,) 'Do you wish fix pressure at sea level (y)',
742	& ' or not (n) ?'
743	read(*,fmt='(a)') yes
744	end do
745
746	if (yes.eq.'y') then
747	write(,) 'Value?'
748	read(,,iostat=ierr) psea
749	DO i=1,iip1
750	DO j=1,jjp1
751	ps(i,j)=psea
752
753	ENDDO
754	ENDDO
755	write(,) 'psea=',psea
756	else
757	write(,) 'no'
758	end if
759
760
761	c emis : change surface emissivity (added by RW)
762	c ----------------------------------------------
763	else if (trim(modif).eq.'emis') then
764
765	print*,'new value?'
766	read(,) emisread
767
768	do i=1,ngridmx
769	emis(i)=emisread
770	enddo
771
772	c qname : change tracer name
773	c --------------------------
774	else if (trim(modif).eq.'qname') then
775	yes='y'
776	do while (yes.eq.'y')
777	write(,) 'Which tracer name do you want to change ?'
778	do iq=1,nqtot
779	write(*,'(i3,a3,a20)')iq,' : ',trim(tname(iq))
780	enddo
781	write(*,'(a35,i3)')
782	& '(enter tracer number; between 1 and ',nqtot
783	write(,)' or any other value to quit this option)'
784	read(,) iq
785	if ((iq.ge.1).and.(iq.le.nqtot)) then
786	write(,)'Change tracer name ',trim(tname(iq)),' to ?'
787	read(,) txt
788	tname(iq)=txt
789	write(,)'Do you want to change another tracer name (y/n)?'
790	read(*,'(a)') yes
791	else
792	! inapropiate value of iq; quit this option
793	yes='n'
794	endif ! of if ((iq.ge.1).and.(iq.le.nqtot))
795	enddo ! of do while (yes.ne.'y')
796
797	c q=0 : set tracers to zero
798	c -------------------------
799	else if (trim(modif).eq.'q=0') then
800	c mise a 0 des q (traceurs)
801	write(,) 'Tracers set to 0 (1.E-30 in fact)'
802	DO iq =1, nqtot
803	DO l=1,llm
804	DO j=1,jjp1
805	DO i=1,iip1
806	q(i,j,l,iq)=1.e-30
807	ENDDO
808	ENDDO
809	ENDDO
810	ENDDO
811
812	c set surface tracers to zero
813	DO iq =1, nqtot
814	DO ig=1,ngridmx
815	qsurf(ig,iq)=0.
816	ENDDO
817	ENDDO
818
819	c q=x : initialise tracer manually
820	c --------------------------------
821	else if (trim(modif).eq.'q=x') then
822	write(,) 'Which tracer do you want to modify ?'
823	do iq=1,nqtot
824	write(,)iq,' : ',trim(tname(iq))
825	enddo
826	write(,) '(choose between 1 and ',nqtot,')'
827	read(,) iq
828	write(,)'mixing ratio of tracer ',trim(tname(iq)),
829	& ' ? (kg/kg)'
830	read(,) val
831	DO l=1,llm
832	DO j=1,jjp1
833	DO i=1,iip1
834	q(i,j,l,iq)=val
835	ENDDO
836	ENDDO
837	ENDDO
838	write(,) 'SURFACE value of tracer ',trim(tname(iq)),
839	& ' ? (kg/m2)'
840	read(,) val
841	DO ig=1,ngridmx
842	qsurf(ig,iq)=val
843	ENDDO
844
845	c t=profile : initialize temperature with a given profile
846	c -------------------------------------------------------
847	else if (trim(modif) .eq. 't=profile') then
848	write(,) 'Temperature profile from ASCII file'
849	write(,) "'profile.in' e.g. 1D output"
850	write(,) "(one value per line in file; starting with"
851	write(,) "surface value, the 1st atmospheric layer"
852	write(,) "followed by 2nd, etc. up to top of atmosphere)"
853	txt="profile.in"
854	open(33,file=trim(txt),status='old',form='formatted',
855	& iostat=ierr)
856	if (ierr.eq.0) then
857	! OK, found file 'profile_...', load the profile
858	do l=1,llm+1
859	read(33,*,iostat=ierr) profile(l)
860	write(,) profile(l)
861	if (ierr.ne.0) then ! something went wrong
862	exit ! quit loop
863	endif
864	enddo
865	if (ierr.eq.0) then
866	tsurf(1:ngridmx)=profile(1)
867	tsoil(1:ngridmx,1:nsoilmx)=profile(1)
868	do l=1,llm
869	Tset(1:iip1,1:jjp1,l)=profile(l+1)
870	flagtset=.true.
871	enddo
872	ucov(1:iip1,1:jjp1,1:llm)=0.
873	vcov(1:iip1,1:jjm,1:llm)=0.
874	q2(1:ngridmx,1:llm+1)=0.
875	else
876	write(,)'problem reading file ',trim(txt),' !'
877	write(,)'No modifications to temperature'
878	endif
879	else
880	write(,)'Could not find file ',trim(txt),' !'
881	write(,)'No modifications to temperature'
882	endif
883
884	c q=profile : initialize tracer with a given profile
885	c --------------------------------------------------
886	else if (trim(modif) .eq. 'q=profile') then
887	write(,) 'Tracer profile will be sought in ASCII file'
888	write(,) "'profile_tracer' where 'tracer' is tracer name"
889	write(,) "(one value per line in file; starting with"
890	write(,) "surface value, the 1st atmospheric layer"
891	write(,) "followed by 2nd, etc. up to top of atmosphere)"
892	write(,) 'Which tracer do you want to set?'
893	do iq=1,nqtot
894	write(,)iq,' : ',trim(tname(iq))
895	enddo
896	write(,) '(choose between 1 and ',nqtot,')'
897	read(,) iq
898	if ((iq.lt.1).or.(iq.gt.nqtot)) then
899	! wrong value for iq, go back to menu
900	write(,) "wrong input value:",iq
901	cycle
902	endif
903	! look for input file 'profile_tracer'
904	txt="profile_"//trim(tname(iq))
905	open(41,file=trim(txt),status='old',form='formatted',
906	& iostat=ierr)
907	if (ierr.eq.0) then
908	! OK, found file 'profile_...', load the profile
909	do l=1,llm+1
910	read(41,*,iostat=ierr) profile(l)
911	if (ierr.ne.0) then ! something went wrong
912	exit ! quit loop
913	endif
914	enddo
915	if (ierr.eq.0) then
916	! initialize tracer values
917	qsurf(:,iq)=profile(1)
918	do l=1,llm
919	q(:,:,l,iq)=profile(l+1)
920	enddo
921	write(,)'OK, tracer ',trim(tname(iq)),' initialized '
922	& ,'using values from file ',trim(txt)
923	else
924	write(,)'problem reading file ',trim(txt),' !'
925	write(,)'No modifications to tracer ',trim(tname(iq))
926	endif
927	else
928	write(,)'Could not find file ',trim(txt),' !'
929	write(,)'No modifications to tracer ',trim(tname(iq))
930	endif
931
932
933	c isotherm : Isothermal temperatures and no winds
934	c -----------------------------------------------
935	else if (trim(modif) .eq. 'isotherm') then
936
937	write(,)'Isothermal temperature of the atmosphere,
938	& surface and subsurface'
939	write(,) 'Value of this temperature ? :'
940	203 read(,,iostat=ierr) Tiso
941	if(ierr.ne.0) goto 203
942
943	tsurf(1:ngridmx)=Tiso
944
945	tsoil(1:ngridmx,1:nsoilmx)=Tiso
946
947	Tset(1:iip1,1:jjp1,1:llm)=Tiso
948	flagtset=.true.
949
950	t(1:iip1,1:jjp1,1:llm)=Tiso
951	!! otherwise hydrost. integrations below
952	!! use the wrong temperature
953	!! -- NB: Tset might be useless
954
955	ucov(1:iip1,1:jjp1,1:llm)=0
956	vcov(1:iip1,1:jjm,1:llm)=0
957	q2(1:ngridmx,1:llm+1)=0
958
959	c radequi : Radiative equilibrium profile of temperatures and no winds
960	c --------------------------------------------------------------------
961	else if (trim(modif) .eq. 'radequi') then
962
963	write(,)'radiative equilibrium temperature profile'
964
965	do ig=1, ngridmx
966	teque= 335.0-60.0sin(latfi(ig))sin(latfi(ig))-
967	& 10.0cos(latfi(ig))cos(latfi(ig))
968	tsurf(ig) = MAX(220.0,teque)
969	end do
970	do l=2,nsoilmx
971	do ig=1, ngridmx
972	tsoil(ig,l) = tsurf(ig)
973	end do
974	end do
975	DO j=1,jjp1
976	DO i=1,iim
977	Do l=1,llm
978	teque=335.-60.0sin(rlatu(j))sin(rlatu(j))-
979	& 10.0cos(rlatu(j))cos(rlatu(j))
980	Tset(i,j,l)=MAX(220.0,teque)
981	end do
982	end do
983	end do
984	flagtset=.true.
985	ucov(1:iip1,1:jjp1,1:llm)=0
986	vcov(1:iip1,1:jjm,1:llm)=0
987	q2(1:ngridmx,1:llm+1)=0
988
989	! therm_ini_s: (re)-set soil thermal inertia to reference surface values
990	! ----------------------------------------------------------------------
991
992	else if (trim(modif) .eq. 'therm_ini_s') then
993	! write(,)"surfithfi(1):",surfithfi(1)
994	do isoil=1,nsoilmx
995	inertiedat(1:ngridmx,isoil)=surfithfi(1:ngridmx)
996	enddo
997	write(,)'OK: Soil thermal inertia has been reset to referenc
998	&e surface values'
999	! write(,)"inertiedat(1,1):",inertiedat(1,1)
1000	ithfi(:,:)=inertiedat(:,:)
1001	! recast ithfi() onto ith()
1002	call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
1003	! Check:
1004	! do i=1,iip1
1005	! do j=1,jjp1
1006	! do isoil=1,nsoilmx
1007	! write(77,*) i,j,isoil," ",ith(i,j,isoil)
1008	! enddo
1009	! enddo
1010	! enddo
1011
1012
1013
1014	c slab_ocean_initialisation
1015	c ------------------------------------------------
1016	else if (trim(modif) .eq. 'slab_ocean_0') then
1017
1018	write(,) "No ocean yet on Titan ! Can't use this option"
1019	stop
1020
1021	else
1022	write(,) ' Unknown (misspelled?) option!!!'
1023	end if ! of if (trim(modif) .eq. '...') elseif ...
1024
1025
1026
1027	enddo ! of do ! infinite loop on liste of changes
1028
1029	999 continue
1030
1031
1032	c================================================================
1033	c Initialisation for methane surface tank (added by JVO 2017)
1034	c================================================================
1035	DO ig=1, ngridmx
1036	tankCH4(ig)=2.0
1037	ENDDO
1038
1039	c=======================================================================
1040	c Correct pressure on the new grid (menu 0)
1041	c=======================================================================
1042
1043
1044	if ((choix_1.eq.0).and.(.not.flagps0)) then
1045	r = 1000.*8.31/mugaz
1046
1047	do j=1,jjp1
1048	do i=1,iip1
1049	ps(i,j) = ps(i,j) *
1050	. exp((phisold_newgrid(i,j)-phis(i,j)) /
1051	. (t(i,j,1) * r))
1052	end do
1053	end do
1054
1055	c periodicite de ps en longitude
1056	do j=1,jjp1
1057	ps(1,j) = ps(iip1,j)
1058	end do
1059	end if
1060
1061
1062	c=======================================================================
1063	c=======================================================================
1064
1065	c=======================================================================
1066	c Initialisation de la physique / ecriture de newstartfi :
1067	c=======================================================================
1068
1069
1070	CALL inifilr
1071	CALL pression(ip1jmp1, ap, bp, ps, p3d)
1072
1073
1074	c=========================================================================
1075	c Calcul de la dimension verticale pour la chimie - JVO 2017
1076	c start_archive seulement, la grille verticale pouvant avoir ete modifiee
1077	c==========================================================================
1078
1079	IF (choix_1.eq.0) THEN
1080
1081	! Calculate the # of upper chemistry layers with the "new" pressure grid
1082	! For this we use Vervack profile for upper atmosphere with dz=10km
1083
1084	CALL gr_kim_vervack
1085
1086	WRITE(,)
1087	WRITE(,) " With the compiled vertical grid we found :"
1088	WRITE(,) " Number of upper chemistry layers =", nlaykim_up
1089
1090	IF (callchim) THEN
1091
1092	IF (.NOT.allocated(ykim_up)) THEN
1093	ALLOCATE(ykim_up(nkim,ngridmx,nlaykim_up))
1094	ENDIF
1095
1096	! Regridind of upper chemistry fields and uppermost GCM
1097	! chemistry tracers if needed
1098
1099	IF (nlaykimold.ne.nlaykim_up) THEN
1100
1101	WRITE(,) " Warning, nlaykimold=", nlaykimold
1102	WRITE(,) ' which implies that a vertical regrid on upper
1103	&chemistry fields will be performed. If the ceiling is lowered I
1104	&will also adjust GCM uppermost layers chem. tracers.'
1105	WRITE(,)
1106
1107	CALL vert_regrid_kim(nqtot,q)
1108
1109	ELSE
1110	ykim_up(:,:,:) = ykim_up_oldv(:,:,:) ! No vertical recast
1111	ENDIF
1112
1113	ENDIF
1114
1115	endif ! of if (choix_1.eq.0)
1116
1117
1118	c-----------------------------------------------------------------------
1119	c Initialisation pks:
1120	c-----------------------------------------------------------------------
1121
1122	CALL exner_hyb(ip1jmp1, ps, p3d, pks, pk, pkf)
1123	! Calcul de la temperature potentielle teta
1124
1125	if (flagtset) then
1126	DO l=1,llm
1127	DO j=1,jjp1
1128	DO i=1,iim
1129	teta(i,j,l) = Tset(i,j,l) * cpp/pk(i,j,l)
1130	ENDDO
1131	teta (iip1,j,l)= teta (1,j,l)
1132	ENDDO
1133	ENDDO
1134	else if (choix_1.eq.0) then
1135	DO l=1,llm
1136	DO j=1,jjp1
1137	DO i=1,iim
1138	teta(i,j,l) = t(i,j,l) * cpp/pk(i,j,l)
1139	ENDDO
1140	teta (iip1,j,l)= teta (1,j,l)
1141	ENDDO
1142	ENDDO
1143	endif
1144
1145	C Calcul intermediaire
1146	c
1147	if (choix_1.eq.0) then
1148	CALL massdair( p3d, masse )
1149	c
1150	print *,' ALPHAX ',alphax
1151
1152	DO l = 1, llm
1153	DO i = 1, iim
1154	xppn(i) = aire( i, 1 ) * masse( i , 1 , l )
1155	xpps(i) = aire( i,jjp1 ) * masse( i , jjp1 , l )
1156	ENDDO
1157	xpn = SUM(xppn)/apoln
1158	xps = SUM(xpps)/apols
1159	DO i = 1, iip1
1160	masse( i , 1 , l ) = xpn
1161	masse( i , jjp1 , l ) = xps
1162	ENDDO
1163	ENDDO
1164	endif
1165	phis(iip1,:) = phis(1,:)
1166
1167	itau=0
1168	if (choix_1.eq.0) then
1169	day_ini=int(date)
1170	endif
1171
1172	c
1173	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
1174
1175	CALL caldyn0( itau,ucov,vcov,teta,ps,masse,pk,phis ,
1176	* phi,w, pbaru,pbarv,day_ini+time )
1177
1178
1179	CALL dynredem0("restart.nc",day_ini,phis)
1180	CALL dynredem1("restart.nc",0.0,vcov,ucov,teta,q,masse,ps)
1181	C
1182	C Ecriture etat initial physique
1183	C
1184
1185	call physdem0("restartfi.nc",lonfi,latfi,nsoilmx,ngridmx,llm,
1186	& nqtot,dtphys,real(day_ini),0.0,
1187	& airefi,albfi,ithfi,zmea,zstd,zsig,zgam,zthe)
1188	call physdem1("restartfi.nc",nsoilmx,ngridmx,llm,nqtot,
1189	& dtphys,real(day_ini),
1190	& tsurf,tsoil,emis,q2,qsurf,tankCH4)
1191
1192	c=======================================================================
1193	c Formats
1194	c=======================================================================
1195
1196	1 FORMAT(//10x,'la valeur de im =',i4,2x,'lue sur le fichier de dema
1197	*rrage est differente de la valeur parametree iim =',i4//)
1198	2 FORMAT(//10x,'la valeur de jm =',i4,2x,'lue sur le fichier de dema
1199	*rrage est differente de la valeur parametree jjm =',i4//)
1200	3 FORMAT(//10x,'la valeur de lllm =',i4,2x,'lue sur le fichier demar
1201	*rage est differente de la valeur parametree llm =',i4//)
1202
1203	write(,) "newstart: All is well that ends well."
1204
1205	end
1206

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

Download in other formats:

Original Format