Context Navigation

newstart.F @ 787

Last change on this file since 787 was 787, checked in by aslmd, 12 years ago

LMDZ.GENERIC. (Sorry for long text but this is a quite major commit)

Paving the path for parallel computations. And moving towards a more flexible code.

Automatic allocation is used within all routines in phystd. No further mention to ngridmx and nqmx.

ngridmx and nqmx are still used in LMDZ.GENERIC in the dyn3d part
if the LMDZ4/LMDZ5 dynamical core is used, there is no more fixed dimensions ngridmx and nqmx --> a fully flexible parallel implementation is now possible (e.g. no need to recompile when changing numbers of processors)

The important stuff :

Compilation checked with ifort. OK with and without debug mode. No errors. Checked for: gcm, newstart, rcm1d, kcm1d
RUN GCM: Running an Earth test case. Comparison with previous revision --> debug mode : perfect match. bit by bit (diff command). checked with plots --> O1 mode : close match (checked with plots) --> O2 mode : sometimes up to 0.5 K departure.... BUT in this new version O2 and O1 are quite close while in previous version O1 and O2 differed by about, well, typically 0.5 K (pictures available on request)
RUN NEWSTART : perfect match (bit-by-bit) in either debug or normal mode.
RUN RCM1D : perfect match in normal mode.
RUN KCM1D : not tested (I don't know what is the use of kcm1d)

List of main changes :

Additional arguments to some subroutines (ngrid and nq)
F77 include strategy is obsolete and replaced by F90 module strategy In this new strategy arrays are allocatable and allocated once at first use This has to be done for all common featuring arrays defined with ngridmx or nqmx

surfdat.h >> surfdat_h.F90
tracer.h >> tracer_h.F90
comsaison.h >> comsaison_h.F90
comgeomfi.h >> comgeomfi_h.F90
comsoil.h >> comsoil_h.F90
comdiurn.h >> comdiurn_h.F90
fisice.h >> DELETED. was not used. probably a fossil.
watercap.h >> DELETED. variable put in surfdat_h.F90

F77 'save' strategy is obsolete and replaced by F90 'allocatable save' strategy (see previous point and e.g. new version of physiq.F90)
Suppressing any mention to advtrac.h which is a common in the dynamics and needs nqmx This was easily solved by adding an argument with tracer names, coming from the dynamics This is probably not a definitive solution, ... but this allows for generic physics to work easily with either LMDZ.GENERIC or LMDZ dynamical cores
Removing consistency tests between nq and nqmx ; and ngrid and ngridmx. No use now!
Adaptation of rcm1d, kcm1d, newstart given above-mentioned changes

A note on phyetat0 and soil_setting:

Now written so that a slice of horizontal size 'ngrid' starting at grid point 'cursor' is read in startfi.nc 'cursor' is defined in dimphys.h and initialized by inifis (or in newstart) this is useful for parallel computations. default behavior is the usual one : sequential runs, cursor is 1, size ngrid is the whole global domain

A note on an additional change :

nueffrad is now an argument to callcorrk as is the case for reffrad both are saved in physiq this is for consistency and lisibility (previously nueffrad was saved in callcorrk) ... but there is a call to a function which modifies nueffrad in physiq ... previously this was not modifying nueffrad (although it was quite cumbersome to detect this) ... to be conservative I kept this behaviour and highlighted it with an array nueffrad_dummy ... I added a comment because someone might want to change this

File size: 48.3 KB

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

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

Download in other formats:

Original Format