Context Navigation

newstart.F @ 3624

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

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

Download in other formats:

Original Format