Context Navigation

newstart.F @ 1150

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

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

Download in other formats:

Original Format