Context Navigation

newstart.F @ 1644

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

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

Download in other formats:

Original Format