Context Navigation

newstart.F @ 1230

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

Generic model:
Major cleanup, in order to ease the use of LMDZ.GENERIC with (parallel) dynamics
in LMDZ.COMMON: (NB: this will break LMDZ.UNIVERSAL, which should be thrashed
in the near future)

Updated makegcm_* scripts (and makdim) and added the "-full" (to enforce full recomputation of the model) option
In dyn3d: converted control.h to module control_mod.F90 and converted iniadvtrac.F to module infotrac.F90
Added module mod_const_mpi.F90 in dyn3d (not used in serial mode)
Rearanged input/outputs routines everywhere to handle serial/MPI cases. physdem.F => phyredem.F90 , phyetat0.F => phyetat0.F90 ; all read/write routines for startfi files are gathered in module iostart.F90
added parallelism related routines init_phys_lmdz.F90, comgeomphy.F90, dimphy.F90, iniphysiq.F90, mod_grid_phy_lmdz.F90, mod_phys_lmdz_mpi_data.F90, mod_phys_lmdz_mpi_transfert.F90, mod_phys_lmdz_omp_data.F90, mod_phys_lmdz_omp_transfert.F90, mod_phys_lmdz_para.F90, mod_phys_lmdz_transfert_para.F90 in phymars and mod_const_mpi.F90 in dyn3d (for compliance with parallelism)
added created generic routines 'planetwide_maxval' and 'planetwide_minval', in module "planetwide_mod", that enable obtaining the max and min of a field over the whole planet. This should be further imroved with computation of means (possibly area weighed), etc.

File size: 53.6 KB

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

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

Download in other formats:

Original Format