Context Navigation

newstart.F @ 1538

Last change on this file since 1538 was 1538, checked in by emillour, 9 years ago

Generic GCM:

Made nsoilmx be no longer a "parameter" and thus added the possibility to define the number of subsurface layers nsoilmx, along with first layer thickness "lay1_soil" and companion coefficient "alpha_soil", in callphys.def at run time. As before (when these were hard-coded), these are such that the depth of soil mid-layers are: mlayer(k)=lay1_soil*alpha_soil(k-1/2), for k=0,nsoil-1

File size: 58.0 KB

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

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

Download in other formats:

Original Format