Context Navigation

newstart.F @ 1328

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

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

Download in other formats:

Original Format