Context Navigation

newstart.F @ 1252

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

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

Download in other formats:

Original Format