Context Navigation

newstart.F @ 988

Last change on this file since 988 was 988, checked in by emillour, 12 years ago

Generic GCM:

Moved "newstart" (and related "lect_start_archive.F") to phystd directory
Adapted makegcm_* scripts to enable compiling main prog from physics
Added in newstart the possibility to not read in any surface.nc file (when loading a start_archive) with keyword "none" (instead of surface file name)
Some general cleanup:
- in bibio: removed unused lmdstd.h readstd.F writestd.F mywrite.F

readcoord.F scatter.F gather.F ini36.F from36.F to36.F
lnblnk.F (F90 len_trim() should be used instead)

in dyn3d: removed unused inigrads.F wrgrads.F gradsdef.h

xvik.F (specific to current Mars)

File size: 50.9 KB

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

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

Download in other formats:

Original Format