Context Navigation

newstart.F @ 2613

Last change on this file since 2613 was 2354, checked in by emillour, 5 years ago
Generic GCM: Major cleanup: remove obsolete compilation scripts (makegcm*) and old dynamical core, as it is obsolete with respect to the one provide in LMDZ.COMMON. 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	call initracer(ngridmx,nqtot,tname)
288
289	! Initialize dimphy module (klon,klev,..)
290	call init_dimphy(ngridmx,llm)
291	! Allocate saved arrays (as in firstcall of physiq)
292	call phys_state_var_init(nqtot)
293
294	c-----------------------------------------------------------------------
295	c Lecture du tableau des parametres du run (pour la dynamique)
296	c-----------------------------------------------------------------------
297	if (choix_1.eq.0) then
298
299	write(,) 'reading tab_cntrl START_ARCHIVE'
300	c
301	ierr = NF_INQ_VARID (nid, "controle", nvarid)
302	#ifdef NC_DOUBLE
303	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tab_cntrl)
304	#else
305	ierr = NF_GET_VAR_REAL(nid, nvarid, tab_cntrl)
306	#endif
307	c
308	else if (choix_1.eq.1) then
309
310	write(,) 'reading tab_cntrl START'
311	c
312	ierr = NF_INQ_VARID (nid_dyn, "controle", nvarid)
313	#ifdef NC_DOUBLE
314	ierr = NF_GET_VAR_DOUBLE(nid_dyn, nvarid, tab_cntrl)
315	#else
316	ierr = NF_GET_VAR_REAL(nid_dyn, nvarid, tab_cntrl)
317	#endif
318	c
319	write(,) 'reading tab_cntrl STARTFI'
320	c
321	ierr = NF_INQ_VARID (nid_fi, "controle", nvarid)
322	#ifdef NC_DOUBLE
323	ierr = NF_GET_VAR_DOUBLE(nid_fi, nvarid, tab_cntrl_bis)
324	#else
325	ierr = NF_GET_VAR_REAL(nid_fi, nvarid, tab_cntrl_bis)
326	#endif
327	c
328	do i=1,50
329	tab_cntrl(i+50)=tab_cntrl_bis(i)
330	enddo
331	write(,) 'printing tab_cntrl', tab_cntrl
332	do i=1,100
333	write(,) i,tab_cntrl(i)
334	enddo
335
336	write(,) 'Reading file START'
337	fichnom = 'start.nc'
338	CALL dynetat0(fichnom,vcov,ucov,teta,q,masse,
339	. ps,phis,time)
340
341	CALL iniconst
342	CALL inigeom
343
344	! Initialize the physics
345	CALL iniphysiq(iim,jjm,llm,
346	& (jjm-1)*iim+2,comm_lmdz,
347	& daysec,day_ini,dtphys,
348	& rlatu,rlatv,rlonu,rlonv,
349	& aire,cu,cv,rad,g,r,cpp,
350	& 1)
351
352	! Lmodif set to 0 to disable modifications possibility in phyeta0
353	Lmodif=0
354	write(,) 'Reading file STARTFI'
355	fichnom = 'startfi.nc'
356	CALL phyetat0(.true.,ngridmx,llm,fichnom,tab0,Lmodif,nsoilmx,
357	. nqtot,day_ini,time,
358	. tsurf,tsoil,emis,q2,qsurf, !) ! temporary modif by RDW
359	. cloudfrac,totalfrac,hice,rnat,pctsrf_sic,tslab,tsea_ice,
360	. sea_ice)
361
362	! copy albedo and soil thermal inertia on (local) physics grid
363	do i=1,ngridmx
364	albfi(i) = albedodat(i)
365	do j=1,nsoilmx
366	ithfi(i,j) = inertiedat(i,j)
367	enddo
368	! build a surfithfi(:) using 1st layer of ithfi(:), which might
369	! be needed later on if reinitializing soil thermal inertia
370	surfithfi(i)=ithfi(i,1)
371	enddo
372	! also copy albedo and soil thermal inertia on (local) dynamics grid
373	! so that options below can manipulate either (but must then ensure
374	! to correctly recast things on physics grid)
375	call gr_fi_dyn(1,ngridmx,iip1,jjp1,albfi,alb)
376	call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
377	call gr_fi_dyn(1,ngridmx,iip1,jjp1,surfithfi,surfith)
378
379	endif
380	c-----------------------------------------------------------------------
381	c Initialisation des constantes dynamique
382	c-----------------------------------------------------------------------
383
384	kappa = tab_cntrl(9)
385	etot0 = tab_cntrl(12)
386	ptot0 = tab_cntrl(13)
387	ztot0 = tab_cntrl(14)
388	stot0 = tab_cntrl(15)
389	ang0 = tab_cntrl(16)
390	write(,) "Newstart: kappa,etot0,ptot0,ztot0,stot0,ang0"
391	write(,) kappa,etot0,ptot0,ztot0,stot0,ang0
392
393	! for vertical coordinate
394	preff=tab_cntrl(18) ! reference surface pressure
395	pa=tab_cntrl(17) ! reference pressure at which coord is purely pressure
396	!NB: in start_archive files tab_cntrl(17)=tab_cntrl(18)=0
397	write(,) "Newstart: preff=",preff," pa=",pa
398	yes=' '
399	do while ((yes.ne.'y').and.(yes.ne.'n'))
400	write(,) "Change the values of preff and pa? (y/n)"
401	read(*,fmt='(a)') yes
402	end do
403	if (yes.eq.'y') then
404	write(,)"New value of reference surface pressure preff?"
405	write(,)" (for Mars, typically preff=610)"
406	read(,) preff
407	write(,)"New value of reference pressure pa for purely"
408	write(,)"pressure levels (for hybrid coordinates)?"
409	write(,)" (for Mars, typically pa=20)"
410	read(,) pa
411	endif
412	c-----------------------------------------------------------------------
413	c Lecture du tab_cntrl et initialisation des constantes physiques
414	c - pour start: Lmodif = 0 => pas de modifications possibles
415	c (modif dans le tabfi de readfi + loin)
416	c - pour start_archive: Lmodif = 1 => modifications possibles
417	c-----------------------------------------------------------------------
418	if (choix_1.eq.0) then
419	! tabfi requires that input file be first opened by open_startphy(fichnom)
420	fichnom = 'start_archive.nc'
421	call open_startphy(fichnom)
422	call tabfi (ngridmx,nid,Lmodif,tab0,day_ini,lllm,p_rad,
423	. p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
424	else if (choix_1.eq.1) then
425	fichnom = 'startfi.nc'
426	call open_startphy(fichnom)
427	Lmodif=1 ! Lmodif set to 1 to allow modifications in phyeta0
428	call tabfi (ngridmx,nid_fi,Lmodif,tab0,day_ini,lllm,p_rad,
429	. p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
430	endif
431
432	if (p_omeg .eq. -9999.) then
433	p_omeg = 8.*atan(1.)/p_daysec
434	print*,"new value of omega",p_omeg
435	endif
436
437	rad = p_rad
438	omeg = p_omeg
439	g = p_g
440	cpp = p_cpp
441	mugaz = p_mugaz
442	daysec = p_daysec
443
444	if (p_omeg .eq. -9999.) then
445	p_omeg = 8.*atan(1.)/p_daysec
446	print*,"new value of omega",p_omeg
447	endif
448
449	kappa = 8.3141000./(p_mugazp_cpp) ! added by RDW
450
451	c=======================================================================
452	c INITIALISATIONS DIVERSES
453	c=======================================================================
454	! Load parameters from run.def file
455	CALL defrun_new( 99, .TRUE. )
456	! Initialize dynamics geometry and co. (which, when using start.nc may
457	! have changed because of modifications (tabi, preff, pa) above)
458	CALL iniconst
459	CALL inigeom
460	idum=-1
461	idum=0
462
463	! Initialize the physics for start_archive only
464	IF (choix_1.eq.0) THEN
465	CALL iniphysiq(iim,jjm,llm,
466	& (jjm-1)*iim+2,comm_lmdz,
467	& daysec,day_ini,dtphys,
468	& rlatu,rlatv,rlonu,rlonv,
469	& aire,cu,cv,rad,g,r,cpp,
470	& 1)
471	ENDIF
472
473	c=======================================================================
474	c lecture topographie, albedo, inertie thermique, relief sous-maille
475	c=======================================================================
476
477	if (choix_1.eq.0) then ! for start_archive files,
478	! where an external "surface.nc" file is needed
479
480	c do while((relief(1:3).ne.'mol').AND.(relief(1:3).ne.'pla'))
481	c write(,)
482	c write(,) 'choix du relief (mola,pla)'
483	c write(,) '(Topographie MGS MOLA, plat)'
484	c read(*,fmt='(a3)') relief
485	relief="mola"
486	c enddo
487
488	! First get the correct value of datadir, if not already done:
489	! default 'datadir' is set in "datafile_mod"
490	call getin_p("datadir",datadir)
491	write(,) 'Available surface data files are:'
492	filestring='ls '//trim(datadir)//'/'//
493	& trim(surfdir)//' \| grep .nc'
494	call system(filestring)
495	! but in ye old days these files were in datadir, so scan it as well
496	! for the sake of retro-compatibility
497	filestring='ls '//trim(datadir)//'/'//' \| grep .nc'
498	call system(filestring)
499
500	write(,) ''
501	write(,) 'Please choose the relevant file',
502	& ' (e.g. "surface_mars.nc")'
503	write(,) ' or "none" to not use any (and set planetary'
504	write(,) ' albedo and surface thermal inertia)'
505	read(*,fmt='(a50)') surfacefile
506
507	if (surfacefile.ne."none") then
508
509	CALL datareadnc(relief,surfacefile,phis,alb,surfith,
510	& zmeaS,zstdS,zsigS,zgamS,ztheS)
511	else
512	! specific case when not using a "surface.nc" file
513	phis(:,:)=0
514	zmeaS(:,:)=0
515	zstdS(:,:)=0
516	zsigS(:,:)=0
517	zgamS(:,:)=0
518	ztheS(:,:)=0
519
520	write(,) "Enter value of albedo of the bare ground:"
521	read(,) alb(1,1)
522	alb(:,:)=alb(1,1)
523
524	write(,) "Enter value of thermal inertia of soil:"
525	read(,) surfith(1,1)
526	surfith(:,:)=surfith(1,1)
527
528	endif ! of if (surfacefile.ne."none")
529
530	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
531	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,surfith,surfithfi)
532	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
533	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zmeaS,zmea)
534	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zstdS,zstd)
535	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zsigS,zsig)
536	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zgamS,zgam)
537	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ztheS,zthe)
538
539	endif ! of if (choix_1.eq.0)
540
541
542	c=======================================================================
543	c Lecture des fichiers (start ou start_archive)
544	c=======================================================================
545
546	if (choix_1.eq.0) then
547
548	write(,) 'Reading file START_ARCHIVE'
549	CALL lect_start_archive(ngridmx,llm,
550	& date,tsurf,tsoil,emis,q2,
551	& t,ucov,vcov,ps,teta,phisold_newgrid,q,qsurf,
552	& surfith,nid,
553	& rnat,pctsrf_sic,tslab,tsea_ice,sea_ice,
554	& du_nonoro_gwd,dv_nonoro_gwd,east_gwstress,west_gwstress)
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 watercaps : H20 ice on permanent southern cap
1131	c -------------------------------------------------
1132	else if (trim(modif) .eq. 'watercaps') then
1133	if (igcm_h2o_ice.eq.0) then
1134	write(,) "No water ice tracer! Can't use this option"
1135	stop
1136	endif
1137
1138	DO j=1,jjp1
1139	DO i=1,iim
1140	ig=1+(j-2)*iim +i
1141	if(j.eq.1) ig=1
1142	if(j.eq.jjp1) ig=ngridmx
1143
1144	if (rlatu(j)*180./pi.lt.-80.) then
1145	qsurf(ig,igcm_h2o_ice)=3.4e3
1146	!do isoil=1,nsoilmx
1147	! ith(i,j,isoil) = 18000. ! thermal inertia
1148	!enddo
1149	write(,)' ==> Ice mesh South boundary (deg)= ',
1150	& rlatv(j-1)*180./pi
1151	end if
1152	ENDDO
1153	ENDDO
1154	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
1155
1156
1157	c noacglac : H2O ice across highest terrain
1158	c --------------------------------------------
1159	else if (trim(modif) .eq. 'noacglac') then
1160	if (igcm_h2o_ice.eq.0) then
1161	write(,) "No water ice tracer! Can't use this option"
1162	stop
1163	endif
1164	DO j=1,jjp1
1165	DO i=1,iim
1166	ig=1+(j-2)*iim +i
1167	if(j.eq.1) ig=1
1168	if(j.eq.jjp1) ig=ngridmx
1169
1170	if(phis(i,j).gt.1000.*g)then
1171	alb(i,j) = 0.5 ! snow value
1172	do isoil=1,nsoilmx
1173	ith(i,j,isoil) = 18000. ! thermal inertia
1174	! this leads to rnat set to 'ocean' in physiq.F90
1175	! actually no, because it is soil not surface
1176	enddo
1177	endif
1178	ENDDO
1179	ENDDO
1180	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
1181	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
1182	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
1183
1184
1185
1186	c oborealis : H2O oceans across Vastitas Borealis
1187	c -----------------------------------------------
1188	else if (trim(modif) .eq. 'oborealis') then
1189	if (igcm_h2o_ice.eq.0) then
1190	write(,) "No water ice tracer! Can't use this option"
1191	stop
1192	endif
1193	DO j=1,jjp1
1194	DO i=1,iim
1195	ig=1+(j-2)*iim +i
1196	if(j.eq.1) ig=1
1197	if(j.eq.jjp1) ig=ngridmx
1198
1199	if(phis(i,j).lt.-4000.*g)then
1200	! if( (phis(i,j).lt.-4000.*g)
1201	! & .and. (rlatu(j)*180./pi.lt.0.) )then ! south hemisphere only
1202	! phis(i,j)=-8000.0*g ! proper ocean
1203	phis(i,j)=-4000.0*g ! scanty ocean
1204
1205	alb(i,j) = 0.07 ! oceanic value
1206	do isoil=1,nsoilmx
1207	ith(i,j,isoil) = 18000. ! thermal inertia
1208	! this leads to rnat set to 'ocean' in physiq.F90
1209	! actually no, because it is soil not surface
1210	enddo
1211	endif
1212	ENDDO
1213	ENDDO
1214	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
1215	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
1216	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
1217
1218	c iborealis : H2O ice in Northern plains
1219	c --------------------------------------
1220	else if (trim(modif) .eq. 'iborealis') then
1221	if (igcm_h2o_ice.eq.0) then
1222	write(,) "No water ice tracer! Can't use this option"
1223	stop
1224	endif
1225	DO j=1,jjp1
1226	DO i=1,iim
1227	ig=1+(j-2)*iim +i
1228	if(j.eq.1) ig=1
1229	if(j.eq.jjp1) ig=ngridmx
1230
1231	if(phis(i,j).lt.-4000.*g)then
1232	!qsurf(ig,igcm_h2o_ice)=1.e3
1233	qsurf(ig,igcm_h2o_ice)=241.4 ! to make total 33 kg m^-2
1234	endif
1235	ENDDO
1236	ENDDO
1237	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
1238	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
1239	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
1240
1241
1242	c oceanball : H2O liquid everywhere
1243	c ----------------------------
1244	else if (trim(modif) .eq. 'oceanball') then
1245	if (igcm_h2o_ice.eq.0) then
1246	write(,) "No water ice tracer! Can't use this option"
1247	stop
1248	endif
1249	DO j=1,jjp1
1250	DO i=1,iim
1251	ig=1+(j-2)*iim +i
1252	if(j.eq.1) ig=1
1253	if(j.eq.jjp1) ig=ngridmx
1254
1255	qsurf(ig,igcm_h2o_vap)=0.0 ! for ocean, this is infinite source
1256	qsurf(ig,igcm_h2o_ice)=0.0
1257	alb(i,j) = 0.07 ! ocean value
1258
1259	do isoil=1,nsoilmx
1260	ith(i,j,isoil) = 18000. ! thermal inertia
1261	!ith(i,j,isoil) = 50. ! extremely low for test
1262	! this leads to rnat set to 'ocean' in physiq.F90
1263	enddo
1264
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	c iceball : H2O ice everywhere
1272	c ----------------------------
1273	else if (trim(modif) .eq. 'iceball') then
1274	if (igcm_h2o_ice.eq.0) then
1275	write(,) "No water ice tracer! Can't use this option"
1276	stop
1277	endif
1278	DO j=1,jjp1
1279	DO i=1,iim
1280	ig=1+(j-2)*iim +i
1281	if(j.eq.1) ig=1
1282	if(j.eq.jjp1) ig=ngridmx
1283
1284	qsurf(ig,igcm_h2o_vap)=-50. ! for ocean, this is infinite source
1285	qsurf(ig,igcm_h2o_ice)=50. ! == to 0.05 m of oceanic ice
1286	alb(i,j) = 0.6 ! ice albedo value
1287
1288	do isoil=1,nsoilmx
1289	!ith(i,j,isoil) = 18000. ! thermal inertia
1290	! this leads to rnat set to 'ocean' in physiq.F90
1291	enddo
1292
1293	ENDDO
1294	ENDDO
1295	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
1296	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
1297
1298	c supercontinent : H2O ice everywhere
1299	c ----------------------------
1300	else if (trim(modif) .eq. 'supercontinent') then
1301	write(,) 'Minimum longitude (-180,180)'
1302	read(,) val
1303	write(,) 'Maximum longitude (-180,180)'
1304	read(,) val2
1305	write(,) 'Minimum latitude (-90,90)'
1306	read(,) val3
1307	write(,) 'Maximum latitude (-90,90)'
1308	read(,) val4
1309
1310	do j=1,jjp1
1311	do i=1,iip1
1312	ig=1+(j-2)*iim +i
1313	if(j.eq.1) ig=1
1314	if(j.eq.jjp1) ig=ngridmx
1315
1316	c Supercontinent:
1317	if (((rlatu(j)*180./pi.le.val4).and.
1318	& (rlatu(j)*180./pi.ge.val3).and.
1319	& (rlonv(i)*180./pi.le.val2).and.
1320	& (rlonv(i)*180./pi.ge.val))) then
1321
1322	rnat(ig)=1.
1323	alb(i,j) = 0.3
1324	do isoil=1,nsoilmx
1325	ith(i,j,isoil) = 2000.
1326	enddo
1327	c Ocean:
1328	else
1329	rnat(ig)=0.
1330	alb(i,j) = 0.07
1331	do isoil=1,nsoilmx
1332	ith(i,j,isoil) = 18000.
1333	enddo
1334	end if
1335
1336	enddo
1337	enddo
1338	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
1339	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
1340
1341	c isotherm : Isothermal temperatures and no winds
1342	c -----------------------------------------------
1343	else if (trim(modif) .eq. 'isotherm') then
1344
1345	write(,)'Isothermal temperature of the atmosphere,
1346	& surface and subsurface'
1347	write(,) 'Value of this temperature ? :'
1348	203 read(,,iostat=ierr) Tiso
1349	if(ierr.ne.0) goto 203
1350
1351	tsurf(1:ngridmx)=Tiso
1352
1353	tsoil(1:ngridmx,1:nsoilmx)=Tiso
1354
1355	Tset(1:iip1,1:jjp1,1:llm)=Tiso
1356	flagtset=.true.
1357
1358	t(1:iip1,1:jjp1,1:llm)=Tiso
1359	!! otherwise hydrost. integrations below
1360	!! use the wrong temperature
1361	!! -- NB: Tset might be useless
1362
1363	ucov(1:iip1,1:jjp1,1:llm)=0
1364	vcov(1:iip1,1:jjm,1:llm)=0
1365	q2(1:ngridmx,1:llm+1)=0
1366
1367	c radequi : Radiative equilibrium profile of temperatures and no winds
1368	c --------------------------------------------------------------------
1369	else if (trim(modif) .eq. 'radequi') then
1370
1371	write(,)'radiative equilibrium temperature profile'
1372
1373	do ig=1, ngridmx
1374	teque= 335.0-60.0sin(latfi(ig))sin(latfi(ig))-
1375	& 10.0cos(latfi(ig))cos(latfi(ig))
1376	tsurf(ig) = MAX(220.0,teque)
1377	end do
1378	do l=2,nsoilmx
1379	do ig=1, ngridmx
1380	tsoil(ig,l) = tsurf(ig)
1381	end do
1382	end do
1383	DO j=1,jjp1
1384	DO i=1,iim
1385	Do l=1,llm
1386	teque=335.-60.0sin(rlatu(j))sin(rlatu(j))-
1387	& 10.0cos(rlatu(j))cos(rlatu(j))
1388	Tset(i,j,l)=MAX(220.0,teque)
1389	end do
1390	end do
1391	end do
1392	flagtset=.true.
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 coldstart : T set 1K above CO2 frost point and no winds
1398	c ------------------------------------------------
1399	else if (trim(modif) .eq. 'coldstart') then
1400
1401	write(,)'set temperature of the atmosphere,'
1402	&,'surface and subsurface how many degrees above CO2 frost point?'
1403	204 read(,,iostat=ierr) Tabove
1404	if(ierr.ne.0) goto 204
1405
1406	DO j=1,jjp1
1407	DO i=1,iim
1408	ig=1+(j-2)*iim +i
1409	if(j.eq.1) ig=1
1410	if(j.eq.jjp1) ig=ngridmx
1411	tsurf(ig) = (-3167.8)/(log(.01*ps(i,j))-23.23)+Tabove
1412	END DO
1413	END DO
1414	do l=1,nsoilmx
1415	do ig=1, ngridmx
1416	tsoil(ig,l) = tsurf(ig)
1417	end do
1418	end do
1419	DO j=1,jjp1
1420	DO i=1,iim
1421	Do l=1,llm
1422	pp = aps(l) +bps(l)*ps(i,j)
1423	Tset(i,j,l)=(-3167.8)/(log(.01*pp)-23.23)+Tabove
1424	end do
1425	end do
1426	end do
1427
1428	flagtset=.true.
1429	ucov(1:iip1,1:jjp1,1:llm)=0
1430	vcov(1:iip1,1:jjm,1:llm)=0
1431	q2(1:ngridmx,1:llm+1)=0
1432
1433
1434	c co2ice=0 : remove CO2 polar ice caps'
1435	c ------------------------------------------------
1436	else if (trim(modif) .eq. 'co2ice=0') then
1437	! check that there is indeed a co2_ice tracer ...
1438	if (igcm_co2_ice.ne.0) then
1439	do ig=1,ngridmx
1440	!co2ice(ig)=0
1441	qsurf(ig,igcm_co2_ice)=0
1442	emis(ig)=emis(ngridmx/2)
1443	end do
1444	else
1445	write(,) "Can't remove CO2 ice!! (no co2_ice tracer)"
1446	endif
1447
1448	! therm_ini_s: (re)-set soil thermal inertia to reference surface values
1449	! ----------------------------------------------------------------------
1450
1451	else if (trim(modif) .eq. 'therm_ini_s') then
1452	! write(,)"surfithfi(1):",surfithfi(1)
1453	do isoil=1,nsoilmx
1454	inertiedat(1:ngridmx,isoil)=surfithfi(1:ngridmx)
1455	enddo
1456	write(,)'OK: Soil thermal inertia has been reset to referenc
1457	&e surface values'
1458	! write(,)"inertiedat(1,1):",inertiedat(1,1)
1459	ithfi(:,:)=inertiedat(:,:)
1460	! recast ithfi() onto ith()
1461	call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
1462	! Check:
1463	! do i=1,iip1
1464	! do j=1,jjp1
1465	! do isoil=1,nsoilmx
1466	! write(77,*) i,j,isoil," ",ith(i,j,isoil)
1467	! enddo
1468	! enddo
1469	! enddo
1470
1471
1472
1473	c slab_ocean_initialisation
1474	c ------------------------------------------------
1475	else if (trim(modif) .eq. 'slab_ocean_0') then
1476	write(,)'OK: initialisation of slab ocean'
1477
1478	DO ig=1, ngridmx
1479	rnat(ig)=1.
1480	tslab(ig,1)=0.
1481	tslab(ig,2)=0.
1482	tsea_ice(ig)=0.
1483	sea_ice(ig)=0.
1484	pctsrf_sic(ig)=0.
1485
1486	if(ithfi(ig,1).GT.10000.)then
1487	rnat(ig)=0.
1488	phisfi(ig)=0.
1489	tsea_ice(ig)=273.15-1.8
1490	tslab(ig,1)=tsurf(ig)
1491	tslab(ig,2)=tsurf(ig)!*2./3.+(273.15-1.8)/3.
1492	endif
1493
1494	ENDDO
1495	CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
1496
1497
1498	c chemistry_initialisation
1499	c ------------------------------------------------
1500	else if (trim(modif) .eq. 'chemistry_ini') then
1501	write(,)'OK: initialisation of chemical profiles'
1502
1503
1504	call inichim_newstart(ngridmx, nqtot, q, qsurf, ps,
1505	& 1, 0)
1506
1507	! We want to have the very same value at lon -180 and lon 180
1508	do l = 1,llm
1509	do j = 1,jjp1
1510	do iq = 1,nqtot
1511	q(iip1,j,l,iq) = q(1,j,l,iq)
1512	end do
1513	end do
1514	end do
1515
1516
1517	else
1518	write(,) ' Unknown (misspelled?) option!!!'
1519	end if ! of if (trim(modif) .eq. '...') elseif ...
1520
1521
1522
1523	enddo ! of do ! infinite loop on liste of changes
1524
1525	999 continue
1526
1527
1528	c=======================================================================
1529	c Initialisation for cloud fraction and oceanic ice (added by BC 2010)
1530	c=======================================================================
1531	DO ig=1, ngridmx
1532	totalfrac(ig)=0.5
1533	DO l=1,llm
1534	cloudfrac(ig,l)=0.5
1535	ENDDO
1536	! Ehouarn, march 2012: also add some initialisation for hice
1537	hice(ig)=0.0
1538	ENDDO
1539
1540	c========================================================
1541
1542	! DO ig=1,ngridmx
1543	! IF(tsurf(ig) .LT. 273.16-1.8) THEN
1544	! hice(ig)=(273.16-1.8-tsurf(ig))/(273.16-1.8-240)*1.
1545	! hice(ig)=min(hice(ig),1.0)
1546	! ENDIF
1547	! ENDDO
1548
1549
1550
1551
1552	c=======================================================================
1553	c Correct pressure on the new grid (menu 0)
1554	c=======================================================================
1555
1556
1557	if ((choix_1.eq.0).and.(.not.flagps0)) then
1558	r = 1000.*8.31/mugaz
1559
1560	do j=1,jjp1
1561	do i=1,iip1
1562	ps(i,j) = ps(i,j) *
1563	. exp((phisold_newgrid(i,j)-phis(i,j)) /
1564	. (t(i,j,1) * r))
1565	end do
1566	end do
1567
1568	c periodicite de ps en longitude
1569	do j=1,jjp1
1570	ps(1,j) = ps(iip1,j)
1571	end do
1572	end if
1573
1574
1575	c=======================================================================
1576	c=======================================================================
1577
1578	c=======================================================================
1579	c Initialisation de la physique / ecriture de newstartfi :
1580	c=======================================================================
1581
1582
1583	CALL inifilr
1584	CALL pression(ip1jmp1, ap, bp, ps, p3d)
1585
1586	c-----------------------------------------------------------------------
1587	c Initialisation pks:
1588	c-----------------------------------------------------------------------
1589
1590	CALL exner_hyb(ip1jmp1, ps, p3d, pks, pk, pkf)
1591	! Calcul de la temperature potentielle teta
1592
1593	if (flagtset) then
1594	DO l=1,llm
1595	DO j=1,jjp1
1596	DO i=1,iim
1597	teta(i,j,l) = Tset(i,j,l) * cpp/pk(i,j,l)
1598	ENDDO
1599	teta (iip1,j,l)= teta (1,j,l)
1600	ENDDO
1601	ENDDO
1602	else if (choix_1.eq.0) then
1603	DO l=1,llm
1604	DO j=1,jjp1
1605	DO i=1,iim
1606	teta(i,j,l) = t(i,j,l) * cpp/pk(i,j,l)
1607	ENDDO
1608	teta (iip1,j,l)= teta (1,j,l)
1609	ENDDO
1610	ENDDO
1611	endif
1612
1613	C Calcul intermediaire
1614	c
1615	if (choix_1.eq.0) then
1616	CALL massdair( p3d, masse )
1617	c
1618	print *,' ALPHAX ',alphax
1619
1620	DO l = 1, llm
1621	DO i = 1, iim
1622	xppn(i) = aire( i, 1 ) * masse( i , 1 , l )
1623	xpps(i) = aire( i,jjp1 ) * masse( i , jjp1 , l )
1624	ENDDO
1625	xpn = SUM(xppn)/apoln
1626	xps = SUM(xpps)/apols
1627	DO i = 1, iip1
1628	masse( i , 1 , l ) = xpn
1629	masse( i , jjp1 , l ) = xps
1630	ENDDO
1631	ENDDO
1632	endif
1633	phis(iip1,:) = phis(1,:)
1634
1635	itau=0
1636	if (choix_1.eq.0) then
1637	day_ini=int(date)
1638	endif
1639	c
1640	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
1641
1642	CALL caldyn0( itau,ucov,vcov,teta,ps,masse,pk,phis ,
1643	* phi,w, pbaru,pbarv,day_ini+time )
1644
1645
1646	CALL dynredem0("restart.nc",day_ini,phis)
1647	CALL dynredem1("restart.nc",0.0,vcov,ucov,teta,q,masse,ps)
1648	C
1649	C Ecriture etat initial physique
1650	C
1651
1652	call physdem0("restartfi.nc",lonfi,latfi,nsoilmx,ngridmx,llm,
1653	& nqtot,dtphys,real(day_ini),0.0,
1654	& airefi,albfi,ithfi,zmea,zstd,zsig,zgam,zthe)
1655	call physdem1("restartfi.nc",nsoilmx,ngridmx,llm,nqtot,
1656	& dtphys,real(day_ini),
1657	& tsurf,tsoil,emis,q2,qsurf,
1658	& cloudfrac,totalfrac,hice,
1659	& rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
1660
1661	c=======================================================================
1662	c Formats
1663	c=======================================================================
1664
1665	1 FORMAT(//10x,'la valeur de im =',i4,2x,'lue sur le fichier de dema
1666	*rrage est differente de la valeur parametree iim =',i4//)
1667	2 FORMAT(//10x,'la valeur de jm =',i4,2x,'lue sur le fichier de dema
1668	*rrage est differente de la valeur parametree jjm =',i4//)
1669	3 FORMAT(//10x,'la valeur de lllm =',i4,2x,'lue sur le fichier demar
1670	*rage est differente de la valeur parametree llm =',i4//)
1671
1672	write(,) "newstart: All is well that ends well."
1673
1674	end
1675

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

Download in other formats:

Original Format