Context Navigation

newstart.F @ 2236

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

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

Download in other formats:

Original Format