Context Navigation

newstart.F @ 2884

Last change on this file since 2884 was 2785, checked in by emillour, 3 years ago
Generic PCM: Further seperation between dynamics and physics concerning tracers: Tracer names are extracted from traceur.def via initracer.F90 and no longer transfered from the dynamics to the physics LT+EM
File size: 57.8 KB

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

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

Download in other formats:

Original Format