Context Navigation

newstart.F @ 3981

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

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

Download in other formats:

Original Format