Context Navigation

newstart.F @ 1704

Last change on this file since 1704 was 1669, checked in by emillour, 8 years ago

Generic GCM:
Added possibility to run without a startfi.nc file (mainly usefull for
tests with coupling with dynamico dynamical core):

added flag "startphy_file" flag (.false. if doing an "academic" start on the physics side).
turned phyetat0.F90 into module phyetat0_mod.F90
turned tabfi.F into module tabfi_mod.F90 and added handling of startphy_file==.false. case
extra initializations in physiq_mod for startphy_file==.false. case.

File size: 57.6 KB

Line
1	C======================================================================
2	PROGRAM newstart
3	c=======================================================================
4	c
5	c
6	c Auteur: Christophe Hourdin/Francois Forget/Yann Wanherdrick
7	c ------
8	c 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(,) 'q=profile : specify a profile for a tracer'
591	! write(,) 'ini_q : tracers initialisation for chemistry, water an
592	! $d ice '
593	! write(,) 'ini_q-H2O : tracers initialisation for chemistry and
594	! $ice '
595	! write(,) 'ini_q-iceH2O : tracers initialisation for chemistry on
596	! $ly '
597	write(,) 'noglacier : Remove tropical H2O ice if \|lat\|<45'
598	write(,) 'watercapn : H20 ice on permanent N polar cap '
599	write(,) 'watercaps : H20 ice on permanent S polar cap '
600	write(,) 'noacglac : H2O ice across Noachis Terra'
601	write(,) 'oborealis : H2O ice across Vastitas Borealis'
602	write(,) 'iceball : Thick ice layer all over surface'
603	write(,) 'supercontinent: Create a continent of given Ab and TI'
604	write(,) 'wetstart : start with a wet atmosphere'
605	write(,) 'isotherm : Isothermal Temperatures, wind set to zero'
606	write(,) 'radequi : Earth-like radiative equilibrium temperature
607	$ profile (lat-alt) and winds set to zero'
608	write(,) 'coldstart : Start X K above the CO2 frost point and
609	$set wind to zero (assumes 100% CO2)'
610	write(,) 'co2ice=0 : remove CO2 polar cap'
611	write(,) 'ptot : change total pressure'
612	write(,) 'emis : change surface emissivity'
613	write(,) 'therm_ini_s : Set soil thermal inertia to reference sur
614	&face values'
615	write(,) 'slab_ocean_0 : initialisation of slab
616	$ocean variables'
617
618	write(,)
619	write(,) 'Change to perform ?'
620	write(,) ' (enter keyword or return to end)'
621	write(,)
622
623	read(*,fmt='(a20)') modif
624	if (modif(1:1) .eq. ' ')then
625	exit ! exit loop on changes
626	endif
627
628	write(,)
629	write(,) trim(modif) , ' : '
630
631	c 'flat : no topography ("aquaplanet")'
632	c -------------------------------------
633	if (trim(modif) .eq. 'flat') then
634	c set topo to zero
635	z_reel(1:iip1,1:jjp1)=0
636	phis(1:iip1,1:jjp1)=g*z_reel(1:iip1,1:jjp1)
637	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
638	write(,) 'topography set to zero.'
639	write(,) 'WARNING : the subgrid topography parameters',
640	& ' were not set to zero ! => set calllott to F'
641
642	c Choice of surface pressure
643	yes=' '
644	do while ((yes.ne.'y').and.(yes.ne.'n'))
645	write(,) 'Do you wish to choose homogeneous surface',
646	& 'pressure (y) or let newstart interpolate ',
647	& ' the previous field (n)?'
648	read(*,fmt='(a)') yes
649	end do
650	if (yes.eq.'y') then
651	flagps0=.true.
652	write(,) 'New value for ps (Pa) ?'
653	201 read(,,iostat=ierr) patm
654	if(ierr.ne.0) goto 201
655	write(,) patm
656	if (patm.eq.-9999.) then
657	patm = preff
658	write(,) "we set patm = preff", preff
659	endif
660	write(,)
661	write(,) ' new ps everywhere (Pa) = ', patm
662	write(,)
663	do j=1,jjp1
664	do i=1,iip1
665	ps(i,j)=patm
666	enddo
667	enddo
668	end if
669
670	c 'set_ps_to_preff' : used if changing preff with topo
671	c ----------------------------------------------------
672	else if (trim(modif) .eq. 'set_ps_to_preff') then
673	do j=1,jjp1
674	do i=1,iip1
675	ps(i,j)=preff
676	enddo
677	enddo
678
679	c 'nuketharsis : no tharsis bulge for Early Mars'
680	c ---------------------------------------------
681	else if (trim(modif) .eq. 'nuketharsis') then
682
683	DO j=1,jjp1
684	DO i=1,iim
685	ig=1+(j-2)*iim +i
686	if(j.eq.1) ig=1
687	if(j.eq.jjp1) ig=ngridmx
688
689	fact1=(((rlonv(i)180./pi)+100)*2 +
690	& (rlatu(j)180./pi)2)/65*2
691	fact2=exp( -fact1**2.5 )
692
693	phis(i,j) = phis(i,j) - (phis(i,j)+4000.g)fact2
694
695	! if(phis(i,j).gt.2500.*g)then
696	! if(rlatu(j)*180./pi.gt.-80.)then ! avoid chopping south polar cap
697	! phis(i,j)=2500.*g
698	! endif
699	! endif
700
701	ENDDO
702	ENDDO
703	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
704
705
706	c bilball : uniform albedo, thermal inertia
707	c -----------------------------------------
708	else if (trim(modif) .eq. 'bilball') then
709	write(,) 'constante albedo and iner.therm:'
710	write(,) 'New value for albedo (ex: 0.25) ?'
711	101 read(,,iostat=ierr) alb_bb
712	if(ierr.ne.0) goto 101
713	write(,)
714	write(,) ' uniform albedo (new value):',alb_bb
715	write(,)
716
717	write(,) 'New value for thermal inertia (eg: 247) ?'
718	102 read(,,iostat=ierr) ith_bb
719	if(ierr.ne.0) goto 102
720	write(,) 'uniform thermal inertia (new value):',ith_bb
721	DO j=1,jjp1
722	DO i=1,iip1
723	alb(i,j) = alb_bb ! albedo
724	do isoil=1,nsoilmx
725	ith(i,j,isoil) = ith_bb ! thermal inertia
726	enddo
727	END DO
728	END DO
729	! CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ith,ithfi)
730	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
731	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
732
733	c coldspole : sous-sol de la calotte sud toujours froid
734	c -----------------------------------------------------
735	else if (trim(modif) .eq. 'coldspole') then
736	write(,)'new value for the subsurface temperature',
737	& ' beneath the permanent southern polar cap ? (eg: 141 K)'
738	103 read(,,iostat=ierr) tsud
739	if(ierr.ne.0) goto 103
740	write(,)
741	write(,) ' new value of the subsurface temperature:',tsud
742	c nouvelle temperature sous la calotte permanente
743	do l=2,nsoilmx
744	tsoil(ngridmx,l) = tsud
745	end do
746
747
748	write(,)'new value for the albedo',
749	& 'of the permanent southern polar cap ? (eg: 0.75)'
750	104 read(,,iostat=ierr) albsud
751	if(ierr.ne.0) goto 104
752	write(,)
753
754	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
755	c Option 1: only the albedo of the pole is modified :
756	albfi(ngridmx)=albsud
757	write(,) 'ig=',ngridmx,' albedo perennial cap ',
758	& albfi(ngridmx)
759
760	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
761	c Option 2 A haute resolution : coordonnee de la vrai calotte ~
762	c DO j=1,jjp1
763	c DO i=1,iip1
764	c ig=1+(j-2)*iim +i
765	c if(j.eq.1) ig=1
766	c if(j.eq.jjp1) ig=ngridmx
767	c if ((rlatu(j)*180./pi.lt.-84.).and.
768	c & (rlatu(j)*180./pi.gt.-91.).and.
769	c & (rlonv(i)*180./pi.gt.-91.).and.
770	c & (rlonv(i)*180./pi.lt.0.)) then
771	cc albedo de la calotte permanente fixe a albsud
772	c alb(i,j)=albsud
773	c write(,) 'lat=',rlatu(j)*180./pi,
774	c & ' lon=',rlonv(i)*180./pi
775	cc fin de la condition sur les limites de la calotte permanente
776	c end if
777	c ENDDO
778	c ENDDO
779	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
780
781	c CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
782
783
784	c ptot : Modification of the total pressure: ice + current atmosphere
785	c -------------------------------------------------------------------
786	else if (trim(modif).eq.'ptot') then
787
788	! check if we have a co2_ice surface tracer:
789	if (igcm_co2_ice.eq.0) then
790	write(,) " No surface CO2 ice !"
791	write(,) " only atmospheric pressure will be considered!"
792	endif
793	c calcul de la pression totale glace + atm actuelle
794	patm=0.
795	airetot=0.
796	pcap=0.
797	DO j=1,jjp1
798	DO i=1,iim
799	ig=1+(j-2)*iim +i
800	if(j.eq.1) ig=1
801	if(j.eq.jjp1) ig=ngridmx
802	patm = patm + ps(i,j)*aire(i,j)
803	airetot= airetot + aire(i,j)
804	if (igcm_co2_ice.ne.0) then
805	!pcap = pcap + aire(i,j)co2ice(ig)g
806	pcap = pcap + aire(i,j)qsurf(ig,igcm_co2_ice)g
807	endif
808	ENDDO
809	ENDDO
810	ptoto = pcap + patm
811
812	print*,'Current total pressure at surface (co2 ice + atm) ',
813	& ptoto/airetot
814
815	print*,'new value?'
816	read(,) ptotn
817	ptotn=ptotn*airetot
818	patmn=ptotn-pcap
819	print*,'ptoto,patm,ptotn,patmn'
820	print*,ptoto,patm,ptotn,patmn
821	print*,'Mult. factor for pressure (atm only)', patmn/patm
822	do j=1,jjp1
823	do i=1,iip1
824	ps(i,j)=ps(i,j)*patmn/patm
825	enddo
826	enddo
827
828
829
830	c Correction pour la conservation des traceurs
831	yes=' '
832	do while ((yes.ne.'y').and.(yes.ne.'n'))
833	write(,) 'Do you wish to conserve tracer total mass (y)',
834	& ' or tracer mixing ratio (n) ?'
835	read(*,fmt='(a)') yes
836	end do
837
838	if (yes.eq.'y') then
839	write(,) 'OK : conservation of tracer total mass'
840	DO iq =1, nqtot
841	DO l=1,llm
842	DO j=1,jjp1
843	DO i=1,iip1
844	q(i,j,l,iq)=q(i,j,l,iq)*patm/patmn
845	ENDDO
846	ENDDO
847	ENDDO
848	ENDDO
849	else
850	write(,) 'OK : conservation of tracer mixing ratio'
851	end if
852
853	c Correction pour la pression au niveau de la mer
854	yes=' '
855	do while ((yes.ne.'y').and.(yes.ne.'n'))
856	write(,) 'Do you wish fix pressure at sea level (y)',
857	& ' or not (n) ?'
858	read(*,fmt='(a)') yes
859	end do
860
861	if (yes.eq.'y') then
862	write(,) 'Value?'
863	read(,,iostat=ierr) psea
864	DO i=1,iip1
865	DO j=1,jjp1
866	ps(i,j)=psea
867
868	ENDDO
869	ENDDO
870	write(,) 'psea=',psea
871	else
872	write(,) 'no'
873	end if
874
875
876	c emis : change surface emissivity (added by RW)
877	c ----------------------------------------------
878	else if (trim(modif).eq.'emis') then
879
880	print*,'new value?'
881	read(,) emisread
882
883	do i=1,ngridmx
884	emis(i)=emisread
885	enddo
886
887	c qname : change tracer name
888	c --------------------------
889	else if (trim(modif).eq.'qname') then
890	yes='y'
891	do while (yes.eq.'y')
892	write(,) 'Which tracer name do you want to change ?'
893	do iq=1,nqtot
894	write(*,'(i3,a3,a20)')iq,' : ',trim(tname(iq))
895	enddo
896	write(*,'(a35,i3)')
897	& '(enter tracer number; between 1 and ',nqtot
898	write(,)' or any other value to quit this option)'
899	read(,) iq
900	if ((iq.ge.1).and.(iq.le.nqtot)) then
901	write(,)'Change tracer name ',trim(tname(iq)),' to ?'
902	read(,) txt
903	tname(iq)=txt
904	write(,)'Do you want to change another tracer name (y/n)?'
905	read(*,'(a)') yes
906	else
907	! inapropiate value of iq; quit this option
908	yes='n'
909	endif ! of if ((iq.ge.1).and.(iq.le.nqtot))
910	enddo ! of do while (yes.ne.'y')
911
912	c q=0 : set tracers to zero
913	c -------------------------
914	else if (trim(modif).eq.'q=0') then
915	c mise a 0 des q (traceurs)
916	write(,) 'Tracers set to 0 (1.E-30 in fact)'
917	DO iq =1, nqtot
918	DO l=1,llm
919	DO j=1,jjp1
920	DO i=1,iip1
921	q(i,j,l,iq)=1.e-30
922	ENDDO
923	ENDDO
924	ENDDO
925	ENDDO
926
927	c set surface tracers to zero
928	DO iq =1, nqtot
929	DO ig=1,ngridmx
930	qsurf(ig,iq)=0.
931	ENDDO
932	ENDDO
933
934	c q=x : initialise tracer manually
935	c --------------------------------
936	else if (trim(modif).eq.'q=x') then
937	write(,) 'Which tracer do you want to modify ?'
938	do iq=1,nqtot
939	write(,)iq,' : ',trim(tname(iq))
940	enddo
941	write(,) '(choose between 1 and ',nqtot,')'
942	read(,) iq
943	write(,)'mixing ratio of tracer ',trim(tname(iq)),
944	& ' ? (kg/kg)'
945	read(,) val
946	DO l=1,llm
947	DO j=1,jjp1
948	DO i=1,iip1
949	q(i,j,l,iq)=val
950	ENDDO
951	ENDDO
952	ENDDO
953	write(,) 'SURFACE value of tracer ',trim(tname(iq)),
954	& ' ? (kg/m2)'
955	read(,) val
956	DO ig=1,ngridmx
957	qsurf(ig,iq)=val
958	ENDDO
959
960	c t=profile : initialize temperature with a given profile
961	c -------------------------------------------------------
962	else if (trim(modif) .eq. 't=profile') then
963	write(,) 'Temperature profile from ASCII file'
964	write(,) "'profile.in' e.g. 1D output"
965	write(,) "(one value per line in file; starting with"
966	write(,) "surface value, the 1st atmospheric layer"
967	write(,) "followed by 2nd, etc. up to top of atmosphere)"
968	txt="profile.in"
969	open(33,file=trim(txt),status='old',form='formatted',
970	& iostat=ierr)
971	if (ierr.eq.0) then
972	! OK, found file 'profile_...', load the profile
973	do l=1,llm+1
974	read(33,*,iostat=ierr) profile(l)
975	write(,) profile(l)
976	if (ierr.ne.0) then ! something went wrong
977	exit ! quit loop
978	endif
979	enddo
980	if (ierr.eq.0) then
981	tsurf(1:ngridmx)=profile(1)
982	tsoil(1:ngridmx,1:nsoilmx)=profile(1)
983	do l=1,llm
984	Tset(1:iip1,1:jjp1,l)=profile(l+1)
985	flagtset=.true.
986	enddo
987	ucov(1:iip1,1:jjp1,1:llm)=0.
988	vcov(1:iip1,1:jjm,1:llm)=0.
989	q2(1:ngridmx,1:llm+1)=0.
990	else
991	write(,)'problem reading file ',trim(txt),' !'
992	write(,)'No modifications to temperature'
993	endif
994	else
995	write(,)'Could not find file ',trim(txt),' !'
996	write(,)'No modifications to temperature'
997	endif
998
999	c q=profile : initialize tracer with a given profile
1000	c --------------------------------------------------
1001	else if (trim(modif) .eq. 'q=profile') then
1002	write(,) 'Tracer profile will be sought in ASCII file'
1003	write(,) "'profile_tracer' where 'tracer' is tracer name"
1004	write(,) "(one value per line in file; starting with"
1005	write(,) "surface value, the 1st atmospheric layer"
1006	write(,) "followed by 2nd, etc. up to top of atmosphere)"
1007	write(,) 'Which tracer do you want to set?'
1008	do iq=1,nqtot
1009	write(,)iq,' : ',trim(tname(iq))
1010	enddo
1011	write(,) '(choose between 1 and ',nqtot,')'
1012	read(,) iq
1013	if ((iq.lt.1).or.(iq.gt.nqtot)) then
1014	! wrong value for iq, go back to menu
1015	write(,) "wrong input value:",iq
1016	cycle
1017	endif
1018	! look for input file 'profile_tracer'
1019	txt="profile_"//trim(tname(iq))
1020	open(41,file=trim(txt),status='old',form='formatted',
1021	& iostat=ierr)
1022	if (ierr.eq.0) then
1023	! OK, found file 'profile_...', load the profile
1024	do l=1,llm+1
1025	read(41,*,iostat=ierr) profile(l)
1026	if (ierr.ne.0) then ! something went wrong
1027	exit ! quit loop
1028	endif
1029	enddo
1030	if (ierr.eq.0) then
1031	! initialize tracer values
1032	qsurf(:,iq)=profile(1)
1033	do l=1,llm
1034	q(:,:,l,iq)=profile(l+1)
1035	enddo
1036	write(,)'OK, tracer ',trim(tname(iq)),' initialized '
1037	& ,'using values from file ',trim(txt)
1038	else
1039	write(,)'problem reading file ',trim(txt),' !'
1040	write(,)'No modifications to tracer ',trim(tname(iq))
1041	endif
1042	else
1043	write(,)'Could not find file ',trim(txt),' !'
1044	write(,)'No modifications to tracer ',trim(tname(iq))
1045	endif
1046
1047
1048	c wetstart : wet atmosphere with a north to south gradient
1049	c --------------------------------------------------------
1050	else if (trim(modif) .eq. 'wetstart') then
1051	! check that there is indeed a water vapor tracer
1052	if (igcm_h2o_vap.eq.0) then
1053	write(,) "No water vapour tracer! Can't use this option"
1054	stop
1055	endif
1056	DO l=1,llm
1057	DO j=1,jjp1
1058	DO i=1,iip1
1059	q(i,j,l,igcm_h2o_vap)=150.e-6 * (rlatu(j)+pi/2.) / pi
1060	ENDDO
1061	ENDDO
1062	ENDDO
1063
1064	write(,) 'Water mass mixing ratio at north pole='
1065	* ,q(1,1,1,igcm_h2o_vap)
1066	write(,) '---------------------------south pole='
1067	* ,q(1,jjp1,1,igcm_h2o_vap)
1068
1069	c noglacier : remove tropical water ice (to initialize high res sim)
1070	c --------------------------------------------------
1071	else if (trim(modif) .eq. 'noglacier') then
1072	if (igcm_h2o_ice.eq.0) then
1073	write(,) "No water ice tracer! Can't use this option"
1074	stop
1075	endif
1076	do ig=1,ngridmx
1077	j=(ig-2)/iim +2
1078	if(ig.eq.1) j=1
1079	write(,) 'OK: remove surface ice for \|lat\|<45'
1080	if (abs(rlatu(j)*180./pi).lt.45.) then
1081	qsurf(ig,igcm_h2o_ice)=0.
1082	end if
1083	end do
1084
1085
1086	c watercapn : H20 ice on permanent northern cap
1087	c --------------------------------------------------
1088	else if (trim(modif) .eq. 'watercapn') then
1089	if (igcm_h2o_ice.eq.0) then
1090	write(,) "No water ice tracer! Can't use this option"
1091	stop
1092	endif
1093
1094	DO j=1,jjp1
1095	DO i=1,iim
1096	ig=1+(j-2)*iim +i
1097	if(j.eq.1) ig=1
1098	if(j.eq.jjp1) ig=ngridmx
1099
1100	if (rlatu(j)*180./pi.gt.80.) then
1101	qsurf(ig,igcm_h2o_ice)=3.4e3
1102	!do isoil=1,nsoilmx
1103	! ith(i,j,isoil) = 18000. ! thermal inertia
1104	!enddo
1105	write(,)' ==> Ice mesh North boundary (deg)= ',
1106	& rlatv(j-1)*180./pi
1107	end if
1108	ENDDO
1109	ENDDO
1110	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
1111
1112	c$$$ do ig=1,ngridmx
1113	c$$$ j=(ig-2)/iim +2
1114	c$$$ if(ig.eq.1) j=1
1115	c$$$ if (rlatu(j)*180./pi.gt.80.) then
1116	c$$$
1117	c$$$ qsurf(ig,igcm_h2o_ice)=1.e5
1118	c$$$ qsurf(ig,igcm_h2o_vap)=0.0!1.e5
1119	c$$$
1120	c$$$ write(,) 'ig=',ig,' H2O ice mass (kg/m2)= ',
1121	c$$$ & qsurf(ig,igcm_h2o_ice)
1122	c$$$
1123	c$$$ write(,)' ==> Ice mesh South boundary (deg)= ',
1124	c$$$ & rlatv(j)*180./pi
1125	c$$$ end if
1126	c$$$ enddo
1127
1128	c watercaps : H20 ice on permanent southern cap
1129	c -------------------------------------------------
1130	else if (trim(modif) .eq. 'watercaps') then
1131	if (igcm_h2o_ice.eq.0) then
1132	write(,) "No water ice tracer! Can't use this option"
1133	stop
1134	endif
1135
1136	DO j=1,jjp1
1137	DO i=1,iim
1138	ig=1+(j-2)*iim +i
1139	if(j.eq.1) ig=1
1140	if(j.eq.jjp1) ig=ngridmx
1141
1142	if (rlatu(j)*180./pi.lt.-80.) then
1143	qsurf(ig,igcm_h2o_ice)=3.4e3
1144	!do isoil=1,nsoilmx
1145	! ith(i,j,isoil) = 18000. ! thermal inertia
1146	!enddo
1147	write(,)' ==> Ice mesh South boundary (deg)= ',
1148	& rlatv(j-1)*180./pi
1149	end if
1150	ENDDO
1151	ENDDO
1152	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
1153
1154	c$$$ do ig=1,ngridmx
1155	c$$$ j=(ig-2)/iim +2
1156	c$$$ if(ig.eq.1) j=1
1157	c$$$ if (rlatu(j)*180./pi.lt.-80.) then
1158	c$$$ qsurf(ig,igcm_h2o_ice)=1.e5
1159	c$$$ qsurf(ig,igcm_h2o_vap)=0.0 !1.e5
1160	c$$$
1161	c$$$ write(,) 'ig=',ig,' H2O ice mass (kg/m2)= ',
1162	c$$$ & qsurf(ig,igcm_h2o_ice)
1163	c$$$ write(,)' ==> Ice mesh North boundary (deg)= ',
1164	c$$$ & rlatv(j-1)*180./pi
1165	c$$$ end if
1166	c$$$ enddo
1167
1168
1169	c noacglac : H2O ice across highest terrain
1170	c --------------------------------------------
1171	else if (trim(modif) .eq. 'noacglac') then
1172	if (igcm_h2o_ice.eq.0) then
1173	write(,) "No water ice tracer! Can't use this option"
1174	stop
1175	endif
1176	DO j=1,jjp1
1177	DO i=1,iim
1178	ig=1+(j-2)*iim +i
1179	if(j.eq.1) ig=1
1180	if(j.eq.jjp1) ig=ngridmx
1181
1182	if(phis(i,j).gt.1000.*g)then
1183	alb(i,j) = 0.5 ! snow value
1184	do isoil=1,nsoilmx
1185	ith(i,j,isoil) = 18000. ! thermal inertia
1186	! this leads to rnat set to 'ocean' in physiq.F90
1187	! actually no, because it is soil not surface
1188	enddo
1189	endif
1190	ENDDO
1191	ENDDO
1192	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
1193	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
1194	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
1195
1196
1197
1198	c oborealis : H2O oceans across Vastitas Borealis
1199	c -----------------------------------------------
1200	else if (trim(modif) .eq. 'oborealis') then
1201	if (igcm_h2o_ice.eq.0) then
1202	write(,) "No water ice tracer! Can't use this option"
1203	stop
1204	endif
1205	DO j=1,jjp1
1206	DO i=1,iim
1207	ig=1+(j-2)*iim +i
1208	if(j.eq.1) ig=1
1209	if(j.eq.jjp1) ig=ngridmx
1210
1211	if(phis(i,j).lt.-4000.*g)then
1212	! if( (phis(i,j).lt.-4000.*g)
1213	! & .and. (rlatu(j)*180./pi.lt.0.) )then ! south hemisphere only
1214	! phis(i,j)=-8000.0*g ! proper ocean
1215	phis(i,j)=-4000.0*g ! scanty ocean
1216
1217	alb(i,j) = 0.07 ! oceanic value
1218	do isoil=1,nsoilmx
1219	ith(i,j,isoil) = 18000. ! thermal inertia
1220	! this leads to rnat set to 'ocean' in physiq.F90
1221	! actually no, because it is soil not surface
1222	enddo
1223	endif
1224	ENDDO
1225	ENDDO
1226	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
1227	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
1228	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
1229
1230	c iborealis : H2O ice in Northern plains
1231	c --------------------------------------
1232	else if (trim(modif) .eq. 'iborealis') then
1233	if (igcm_h2o_ice.eq.0) then
1234	write(,) "No water ice tracer! Can't use this option"
1235	stop
1236	endif
1237	DO j=1,jjp1
1238	DO i=1,iim
1239	ig=1+(j-2)*iim +i
1240	if(j.eq.1) ig=1
1241	if(j.eq.jjp1) ig=ngridmx
1242
1243	if(phis(i,j).lt.-4000.*g)then
1244	!qsurf(ig,igcm_h2o_ice)=1.e3
1245	qsurf(ig,igcm_h2o_ice)=241.4 ! to make total 33 kg m^-2
1246	endif
1247	ENDDO
1248	ENDDO
1249	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
1250	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
1251	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
1252
1253
1254	c oceanball : H2O liquid everywhere
1255	c ----------------------------
1256	else if (trim(modif) .eq. 'oceanball') then
1257	if (igcm_h2o_ice.eq.0) then
1258	write(,) "No water ice tracer! Can't use this option"
1259	stop
1260	endif
1261	DO j=1,jjp1
1262	DO i=1,iim
1263	ig=1+(j-2)*iim +i
1264	if(j.eq.1) ig=1
1265	if(j.eq.jjp1) ig=ngridmx
1266
1267	qsurf(ig,igcm_h2o_vap)=0.0 ! for ocean, this is infinite source
1268	qsurf(ig,igcm_h2o_ice)=0.0
1269	alb(i,j) = 0.07 ! ocean value
1270
1271	do isoil=1,nsoilmx
1272	ith(i,j,isoil) = 18000. ! thermal inertia
1273	!ith(i,j,isoil) = 50. ! extremely low for test
1274	! this leads to rnat set to 'ocean' in physiq.F90
1275	enddo
1276
1277	ENDDO
1278	ENDDO
1279	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
1280	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
1281	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
1282
1283	c iceball : H2O ice everywhere
1284	c ----------------------------
1285	else if (trim(modif) .eq. 'iceball') then
1286	if (igcm_h2o_ice.eq.0) then
1287	write(,) "No water ice tracer! Can't use this option"
1288	stop
1289	endif
1290	DO j=1,jjp1
1291	DO i=1,iim
1292	ig=1+(j-2)*iim +i
1293	if(j.eq.1) ig=1
1294	if(j.eq.jjp1) ig=ngridmx
1295
1296	qsurf(ig,igcm_h2o_vap)=-50. ! for ocean, this is infinite source
1297	qsurf(ig,igcm_h2o_ice)=50. ! == to 0.05 m of oceanic ice
1298	alb(i,j) = 0.6 ! ice albedo value
1299
1300	do isoil=1,nsoilmx
1301	!ith(i,j,isoil) = 18000. ! thermal inertia
1302	! this leads to rnat set to 'ocean' in physiq.F90
1303	enddo
1304
1305	ENDDO
1306	ENDDO
1307	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
1308	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
1309
1310	c supercontinent : H2O ice everywhere
1311	c ----------------------------
1312	else if (trim(modif) .eq. 'supercontinent') then
1313	write(,) 'Minimum longitude (-180,180)'
1314	read(,) val
1315	write(,) 'Maximum longitude (-180,180)'
1316	read(,) val2
1317	write(,) 'Minimum latitude (-90,90)'
1318	read(,) val3
1319	write(,) 'Maximum latitude (-90,90)'
1320	read(,) val4
1321
1322	do j=1,jjp1
1323	do i=1,iip1
1324	ig=1+(j-2)*iim +i
1325	if(j.eq.1) ig=1
1326	if(j.eq.jjp1) ig=ngridmx
1327
1328	c Supercontinent:
1329	if (((rlatu(j)*180./pi.le.val4).and.
1330	& (rlatu(j)*180./pi.ge.val3).and.
1331	& (rlonv(i)*180./pi.le.val2).and.
1332	& (rlonv(i)*180./pi.ge.val))) then
1333
1334	rnat(ig)=1.
1335	alb(i,j) = 0.3
1336	do isoil=1,nsoilmx
1337	ith(i,j,isoil) = 2000.
1338	enddo
1339	c Ocean:
1340	else
1341	rnat(ig)=0.
1342	alb(i,j) = 0.07
1343	do isoil=1,nsoilmx
1344	ith(i,j,isoil) = 18000.
1345	enddo
1346	end if
1347
1348	enddo
1349	enddo
1350	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
1351	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
1352
1353	c isotherm : Isothermal temperatures and no winds
1354	c -----------------------------------------------
1355	else if (trim(modif) .eq. 'isotherm') then
1356
1357	write(,)'Isothermal temperature of the atmosphere,
1358	& surface and subsurface'
1359	write(,) 'Value of this temperature ? :'
1360	203 read(,,iostat=ierr) Tiso
1361	if(ierr.ne.0) goto 203
1362
1363	tsurf(1:ngridmx)=Tiso
1364
1365	tsoil(1:ngridmx,1:nsoilmx)=Tiso
1366
1367	Tset(1:iip1,1:jjp1,1:llm)=Tiso
1368	flagtset=.true.
1369
1370	t(1:iip1,1:jjp1,1:llm)=Tiso
1371	!! otherwise hydrost. integrations below
1372	!! use the wrong temperature
1373	!! -- NB: Tset might be useless
1374
1375	ucov(1:iip1,1:jjp1,1:llm)=0
1376	vcov(1:iip1,1:jjm,1:llm)=0
1377	q2(1:ngridmx,1:llm+1)=0
1378
1379	c radequi : Radiative equilibrium profile of temperatures and no winds
1380	c --------------------------------------------------------------------
1381	else if (trim(modif) .eq. 'radequi') then
1382
1383	write(,)'radiative equilibrium temperature profile'
1384
1385	do ig=1, ngridmx
1386	teque= 335.0-60.0sin(latfi(ig))sin(latfi(ig))-
1387	& 10.0cos(latfi(ig))cos(latfi(ig))
1388	tsurf(ig) = MAX(220.0,teque)
1389	end do
1390	do l=2,nsoilmx
1391	do ig=1, ngridmx
1392	tsoil(ig,l) = tsurf(ig)
1393	end do
1394	end do
1395	DO j=1,jjp1
1396	DO i=1,iim
1397	Do l=1,llm
1398	teque=335.-60.0sin(rlatu(j))sin(rlatu(j))-
1399	& 10.0cos(rlatu(j))cos(rlatu(j))
1400	Tset(i,j,l)=MAX(220.0,teque)
1401	end do
1402	end do
1403	end do
1404	flagtset=.true.
1405	ucov(1:iip1,1:jjp1,1:llm)=0
1406	vcov(1:iip1,1:jjm,1:llm)=0
1407	q2(1:ngridmx,1:llm+1)=0
1408
1409	c coldstart : T set 1K above CO2 frost point and no winds
1410	c ------------------------------------------------
1411	else if (trim(modif) .eq. 'coldstart') then
1412
1413	write(,)'set temperature of the atmosphere,'
1414	&,'surface and subsurface how many degrees above CO2 frost point?'
1415	204 read(,,iostat=ierr) Tabove
1416	if(ierr.ne.0) goto 204
1417
1418	DO j=1,jjp1
1419	DO i=1,iim
1420	ig=1+(j-2)*iim +i
1421	if(j.eq.1) ig=1
1422	if(j.eq.jjp1) ig=ngridmx
1423	tsurf(ig) = (-3167.8)/(log(.01*ps(i,j))-23.23)+Tabove
1424	END DO
1425	END DO
1426	do l=1,nsoilmx
1427	do ig=1, ngridmx
1428	tsoil(ig,l) = tsurf(ig)
1429	end do
1430	end do
1431	DO j=1,jjp1
1432	DO i=1,iim
1433	Do l=1,llm
1434	pp = aps(l) +bps(l)*ps(i,j)
1435	Tset(i,j,l)=(-3167.8)/(log(.01*pp)-23.23)+Tabove
1436	end do
1437	end do
1438	end do
1439
1440	flagtset=.true.
1441	ucov(1:iip1,1:jjp1,1:llm)=0
1442	vcov(1:iip1,1:jjm,1:llm)=0
1443	q2(1:ngridmx,1:llm+1)=0
1444
1445
1446	c co2ice=0 : remove CO2 polar ice caps'
1447	c ------------------------------------------------
1448	else if (trim(modif) .eq. 'co2ice=0') then
1449	! check that there is indeed a co2_ice tracer ...
1450	if (igcm_co2_ice.ne.0) then
1451	do ig=1,ngridmx
1452	!co2ice(ig)=0
1453	qsurf(ig,igcm_co2_ice)=0
1454	emis(ig)=emis(ngridmx/2)
1455	end do
1456	else
1457	write(,) "Can't remove CO2 ice!! (no co2_ice tracer)"
1458	endif
1459
1460	! therm_ini_s: (re)-set soil thermal inertia to reference surface values
1461	! ----------------------------------------------------------------------
1462
1463	else if (trim(modif) .eq. 'therm_ini_s') then
1464	! write(,)"surfithfi(1):",surfithfi(1)
1465	do isoil=1,nsoilmx
1466	inertiedat(1:ngridmx,isoil)=surfithfi(1:ngridmx)
1467	enddo
1468	write(,)'OK: Soil thermal inertia has been reset to referenc
1469	&e surface values'
1470	! write(,)"inertiedat(1,1):",inertiedat(1,1)
1471	ithfi(:,:)=inertiedat(:,:)
1472	! recast ithfi() onto ith()
1473	call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
1474	! Check:
1475	! do i=1,iip1
1476	! do j=1,jjp1
1477	! do isoil=1,nsoilmx
1478	! write(77,*) i,j,isoil," ",ith(i,j,isoil)
1479	! enddo
1480	! enddo
1481	! enddo
1482
1483
1484
1485	c slab_ocean_initialisation
1486	c ------------------------------------------------
1487	else if (trim(modif) .eq. 'slab_ocean_0') then
1488	write(,)'OK: initialisation of slab ocean'
1489
1490	DO ig=1, ngridmx
1491	rnat(ig)=1.
1492	tslab(ig,1)=0.
1493	tslab(ig,2)=0.
1494	tsea_ice(ig)=0.
1495	sea_ice(ig)=0.
1496	pctsrf_sic(ig)=0.
1497
1498	if(ithfi(ig,1).GT.10000.)then
1499	rnat(ig)=0.
1500	phisfi(ig)=0.
1501	tsea_ice(ig)=273.15-1.8
1502	tslab(ig,1)=tsurf(ig)
1503	tslab(ig,2)=tsurf(ig)!*2./3.+(273.15-1.8)/3.
1504	endif
1505
1506	ENDDO
1507	CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
1508
1509
1510
1511	else
1512	write(,) ' Unknown (misspelled?) option!!!'
1513	end if ! of if (trim(modif) .eq. '...') elseif ...
1514
1515
1516
1517	enddo ! of do ! infinite loop on liste of changes
1518
1519	999 continue
1520
1521
1522	c=======================================================================
1523	c Initialisation for cloud fraction and oceanic ice (added by BC 2010)
1524	c=======================================================================
1525	DO ig=1, ngridmx
1526	totalfrac(ig)=0.5
1527	DO l=1,llm
1528	cloudfrac(ig,l)=0.5
1529	ENDDO
1530	! Ehouarn, march 2012: also add some initialisation for hice
1531	hice(ig)=0.0
1532	ENDDO
1533
1534	c========================================================
1535
1536	! DO ig=1,ngridmx
1537	! IF(tsurf(ig) .LT. 273.16-1.8) THEN
1538	! hice(ig)=(273.16-1.8-tsurf(ig))/(273.16-1.8-240)*1.
1539	! hice(ig)=min(hice(ig),1.0)
1540	! ENDIF
1541	! ENDDO
1542
1543
1544
1545
1546	c=======================================================================
1547	c Correct pressure on the new grid (menu 0)
1548	c=======================================================================
1549
1550
1551	if ((choix_1.eq.0).and.(.not.flagps0)) then
1552	r = 1000.*8.31/mugaz
1553
1554	do j=1,jjp1
1555	do i=1,iip1
1556	ps(i,j) = ps(i,j) *
1557	. exp((phisold_newgrid(i,j)-phis(i,j)) /
1558	. (t(i,j,1) * r))
1559	end do
1560	end do
1561
1562	c periodicite de ps en longitude
1563	do j=1,jjp1
1564	ps(1,j) = ps(iip1,j)
1565	end do
1566	end if
1567
1568
1569	c=======================================================================
1570	c=======================================================================
1571
1572	c=======================================================================
1573	c Initialisation de la physique / ecriture de newstartfi :
1574	c=======================================================================
1575
1576
1577	CALL inifilr
1578	CALL pression(ip1jmp1, ap, bp, ps, p3d)
1579
1580	c-----------------------------------------------------------------------
1581	c Initialisation pks:
1582	c-----------------------------------------------------------------------
1583
1584	CALL exner_hyb(ip1jmp1, ps, p3d, beta, pks, pk, pkf)
1585	! Calcul de la temperature potentielle teta
1586
1587	if (flagtset) then
1588	DO l=1,llm
1589	DO j=1,jjp1
1590	DO i=1,iim
1591	teta(i,j,l) = Tset(i,j,l) * cpp/pk(i,j,l)
1592	ENDDO
1593	teta (iip1,j,l)= teta (1,j,l)
1594	ENDDO
1595	ENDDO
1596	else if (choix_1.eq.0) then
1597	DO l=1,llm
1598	DO j=1,jjp1
1599	DO i=1,iim
1600	teta(i,j,l) = t(i,j,l) * cpp/pk(i,j,l)
1601	ENDDO
1602	teta (iip1,j,l)= teta (1,j,l)
1603	ENDDO
1604	ENDDO
1605	endif
1606
1607	C Calcul intermediaire
1608	c
1609	if (choix_1.eq.0) then
1610	CALL massdair( p3d, masse )
1611	c
1612	print *,' ALPHAX ',alphax
1613
1614	DO l = 1, llm
1615	DO i = 1, iim
1616	xppn(i) = aire( i, 1 ) * masse( i , 1 , l )
1617	xpps(i) = aire( i,jjp1 ) * masse( i , jjp1 , l )
1618	ENDDO
1619	xpn = SUM(xppn)/apoln
1620	xps = SUM(xpps)/apols
1621	DO i = 1, iip1
1622	masse( i , 1 , l ) = xpn
1623	masse( i , jjp1 , l ) = xps
1624	ENDDO
1625	ENDDO
1626	endif
1627	phis(iip1,:) = phis(1,:)
1628
1629	itau=0
1630	if (choix_1.eq.0) then
1631	day_ini=int(date)
1632	endif
1633	c
1634	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
1635
1636	CALL caldyn0( itau,ucov,vcov,teta,ps,masse,pk,phis ,
1637	* phi,w, pbaru,pbarv,day_ini+time )
1638
1639
1640	CALL dynredem0("restart.nc",day_ini,phis)
1641	CALL dynredem1("restart.nc",0.0,vcov,ucov,teta,q,masse,ps)
1642	C
1643	C Ecriture etat initial physique
1644	C
1645
1646	call physdem0("restartfi.nc",lonfi,latfi,nsoilmx,ngridmx,llm,
1647	& nqtot,dtphys,real(day_ini),0.0,
1648	& airefi,albfi,ithfi,zmea,zstd,zsig,zgam,zthe)
1649	call physdem1("restartfi.nc",nsoilmx,ngridmx,llm,nqtot,
1650	& dtphys,real(day_ini),
1651	& tsurf,tsoil,emis,q2,qsurf,
1652	& cloudfrac,totalfrac,hice,
1653	& rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
1654
1655	c=======================================================================
1656	c Formats
1657	c=======================================================================
1658
1659	1 FORMAT(//10x,'la valeur de im =',i4,2x,'lue sur le fichier de dema
1660	*rrage est differente de la valeur parametree iim =',i4//)
1661	2 FORMAT(//10x,'la valeur de jm =',i4,2x,'lue sur le fichier de dema
1662	*rrage est differente de la valeur parametree jjm =',i4//)
1663	3 FORMAT(//10x,'la valeur de lllm =',i4,2x,'lue sur le fichier demar
1664	*rage est differente de la valeur parametree llm =',i4//)
1665
1666	write(,) "newstart: All is well that ends well."
1667
1668	end
1669

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

Download in other formats:

Original Format