Context Navigation

newstart.F @ 1351

Last change on this file since 1351 was 1330, checked in by aslmd, 10 years ago
LMDZ.GENERIC. a few corrections in edu version (within deftank).
File size: 55.2 KB

Line
1	C======================================================================
2	PROGRAM newstart
3	c=======================================================================
4	c
5	c
6	c Auteur: Christophe Hourdin/Francois Forget/Yann Wanherdrick
7	c ------
8	c Derniere modif : 12/03
9	c
10	c
11	c Objet: Create or modify the initial state for the LMD Mars GCM
12	c ----- (fichiers NetCDF start et startfi)
13	c
14	c
15	c=======================================================================
16
17	use 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	if (p_omeg .eq. -9999.) then
434	p_omeg = 8.*atan(1.)/p_daysec
435	print*,"new value of omega",p_omeg
436	endif
437
438	kappa = 8.3141000./(p_mugazp_cpp) ! added by RDW
439
440	c=======================================================================
441	c INITIALISATIONS DIVERSES
442	c=======================================================================
443	! Load parameters from run.def file
444	CALL defrun_new( 99, .TRUE. )
445	CALL iniconst
446	CALL inigeom
447	idum=-1
448	idum=0
449
450	c Initialisation coordonnees /aires
451	c -------------------------------
452	! Note: rlatu(:) and rlonv(:) are commons defined in "comgeom.h"
453	! rlatu() and rlonv() are given in radians
454	latfi(1)=rlatu(1)
455	lonfi(1)=0.
456	DO j=2,jjm
457	DO i=1,iim
458	latfi((j-2)*iim+1+i)=rlatu(j)
459	lonfi((j-2)*iim+1+i)=rlonv(i)
460	ENDDO
461	ENDDO
462	latfi(ngridmx)=rlatu(jjp1)
463	lonfi(ngridmx)=0.
464
465	! build airefi(), mesh area on physics grid
466	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
467	! Poles are single points on physics grid
468	airefi(1)=airefi(1)*iim
469	airefi(ngridmx)=airefi(ngridmx)*iim
470
471	! also initialize various physics flags/settings which might be needed
472	! (for instance initracer needs to know about some flags, and/or
473	! 'datafile' path may be changed by user)
474	call inifis(ngridmx,llm,nqtot,day_ini,daysec,dtphys,
475	& latfi,lonfi,airefi,rad,g,r,cpp)
476
477	c=======================================================================
478	c lecture topographie, albedo, inertie thermique, relief sous-maille
479	c=======================================================================
480
481	if (choix_1.eq.0) then ! for start_archive files,
482	! where an external "surface.nc" file is needed
483
484	c do while((relief(1:3).ne.'mol').AND.(relief(1:3).ne.'pla'))
485	c write(,)
486	c write(,) 'choix du relief (mola,pla)'
487	c write(,) '(Topographie MGS MOLA, plat)'
488	c read(*,fmt='(a3)') relief
489	relief="mola"
490	c enddo
491
492	! First get the correct value of datadir, if not already done:
493	! default 'datadir' is set in "datafile_mod"
494	call getin_p("datadir",datadir)
495	write(,) 'Available surface data files are:'
496	filestring='ls '//trim(datadir)//' \| grep .nc'
497	call system(filestring)
498
499	write(,) ''
500	write(,) 'Please choose the relevant file',
501	& ' (e.g. "surface_mars.nc")'
502	write(,) ' or "none" to not use any (and set planetary'
503	write(,) ' albedo and surface thermal inertia)'
504	read(*,fmt='(a50)') surfacefile
505
506	if (surfacefile.ne."none") then
507
508	CALL datareadnc(relief,surfacefile,phis,alb,surfith,
509	& zmeaS,zstdS,zsigS,zgamS,ztheS)
510	else
511	! specific case when not using a "surface.nc" file
512	phis(:,:)=0
513	zmeaS(:,:)=0
514	zstdS(:,:)=0
515	zsigS(:,:)=0
516	zgamS(:,:)=0
517	ztheS(:,:)=0
518
519	write(,) "Enter value of albedo of the bare ground:"
520	read(,) alb(1,1)
521	alb(:,:)=alb(1,1)
522
523	write(,) "Enter value of thermal inertia of soil:"
524	read(,) surfith(1,1)
525	surfith(:,:)=surfith(1,1)
526
527	endif ! of if (surfacefile.ne."none")
528
529	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
530	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,surfith,surfithfi)
531	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
532	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zmeaS,zmea)
533	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zstdS,zstd)
534	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zsigS,zsig)
535	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zgamS,zgam)
536	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ztheS,zthe)
537
538	endif ! of if (choix_1.eq.0)
539
540
541	c=======================================================================
542	c Lecture des fichiers (start ou start_archive)
543	c=======================================================================
544
545	if (choix_1.eq.0) then
546
547	write(,) 'Reading file START_ARCHIVE'
548	CALL lect_start_archive(ngridmx,llm,
549	& date,tsurf,tsoil,emis,q2,
550	& t,ucov,vcov,ps,teta,phisold_newgrid,q,qsurf,
551	& surfith,nid,
552	& rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
553	write(,) "OK, read start_archive file"
554	! copy soil thermal inertia
555	ithfi(:,:)=inertiedat(:,:)
556
557	ierr= NF_CLOSE(nid)
558
559	else if (choix_1.eq.1) then
560	!do nothing, start and startfi have already been read
561	else
562	CALL exit(1)
563	endif
564
565	dtvr = daysec/FLOAT(day_step)
566	dtphys = dtvr * FLOAT(iphysiq)
567
568	c=======================================================================
569	c
570	c=======================================================================
571
572	do ! infinite loop on list of changes
573
574	write(,)
575	write(,)
576	write(,) 'List of possible changes :'
577	write(,) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
578	write(,)
579	write(,) 'flat : no topography ("aquaplanet")'
580	write(,) 'nuketharsis : no Tharsis bulge'
581	write(,) 'bilball : uniform albedo and thermal inertia'
582	write(,) 'coldspole : cold subsurface and high albedo at S.pole'
583	write(,) 'qname : change tracer name'
584	write(,) 't=profile : read temperature profile in profile.in'
585	write(,) 'q=0 : ALL tracer =zero'
586	write(,) 'q=x : give a specific uniform value to one tracer'
587	write(,) 'q=profile : specify a profile for a tracer'
588	! write(,) 'ini_q : tracers initialisation for chemistry, water an
589	! $d ice '
590	! write(,) 'ini_q-H2O : tracers initialisation for chemistry and
591	! $ice '
592	! write(,) 'ini_q-iceH2O : tracers initialisation for chemistry on
593	! $ly '
594	write(,) 'noglacier : Remove tropical H2O ice if \|lat\|<45'
595	write(,) 'watercapn : H20 ice on permanent N polar cap '
596	write(,) 'watercaps : H20 ice on permanent S polar cap '
597	write(,) 'noacglac : H2O ice across Noachis Terra'
598	write(,) 'oborealis : H2O ice across Vastitas Borealis'
599	write(,) 'iceball : Thick ice layer all over surface'
600	write(,) 'wetstart : start with a wet atmosphere'
601	write(,) 'isotherm : Isothermal Temperatures, wind set to zero'
602	write(,) 'radequi : Earth-like radiative equilibrium temperature
603	$ profile (lat-alt) and winds set to zero'
604	write(,) 'coldstart : Start X K above the CO2 frost point and
605	$set wind to zero (assumes 100% CO2)'
606	write(,) 'co2ice=0 : remove CO2 polar cap'
607	write(,) 'ptot : change total pressure'
608	write(,) 'emis : change surface emissivity'
609	write(,) 'therm_ini_s : Set soil thermal inertia to reference sur
610	&face values'
611	write(,) 'slab_ocean_0 : initialisation of slab
612	$ocean variables'
613
614	write(,)
615	write(,) 'Change to perform ?'
616	write(,) ' (enter keyword or return to end)'
617	write(,)
618
619	read(*,fmt='(a20)') modif
620	if (modif(1:1) .eq. ' ')then
621	exit ! exit loop on changes
622	endif
623
624	write(,)
625	write(,) trim(modif) , ' : '
626
627	c 'flat : no topography ("aquaplanet")'
628	c -------------------------------------
629	if (trim(modif) .eq. 'flat') then
630	c set topo to zero
631	z_reel(1:iip1,1:jjp1)=0
632	phis(1:iip1,1:jjp1)=g*z_reel(1:iip1,1:jjp1)
633	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
634	write(,) 'topography set to zero.'
635	write(,) 'WARNING : the subgrid topography parameters',
636	& ' were not set to zero ! => set calllott to F'
637
638	c Choice of surface pressure
639	yes=' '
640	do while ((yes.ne.'y').and.(yes.ne.'n'))
641	write(,) 'Do you wish to choose homogeneous surface',
642	& 'pressure (y) or let newstart interpolate ',
643	& ' the previous field (n)?'
644	read(*,fmt='(a)') yes
645	end do
646	if (yes.eq.'y') then
647	flagps0=.true.
648	write(,) 'New value for ps (Pa) ?'
649	201 read(,,iostat=ierr) patm
650	if(ierr.ne.0) goto 201
651	write(,) patm
652	if (patm.eq.-9999.) then
653	patm = preff
654	write(,) "we set patm = preff", preff
655	endif
656	write(,)
657	write(,) ' new ps everywhere (Pa) = ', patm
658	write(,)
659	do j=1,jjp1
660	do i=1,iip1
661	ps(i,j)=patm
662	enddo
663	enddo
664	end if
665
666	c 'nuketharsis : no tharsis bulge for Early Mars'
667	c ---------------------------------------------
668	else if (trim(modif) .eq. 'nuketharsis') then
669
670	DO j=1,jjp1
671	DO i=1,iim
672	ig=1+(j-2)*iim +i
673	if(j.eq.1) ig=1
674	if(j.eq.jjp1) ig=ngridmx
675
676	fact1=(((rlonv(i)180./pi)+100)*2 +
677	& (rlatu(j)180./pi)2)/65*2
678	fact2=exp( -fact1**2.5 )
679
680	phis(i,j) = phis(i,j) - (phis(i,j)+4000.g)fact2
681
682	! if(phis(i,j).gt.2500.*g)then
683	! if(rlatu(j)*180./pi.gt.-80.)then ! avoid chopping south polar cap
684	! phis(i,j)=2500.*g
685	! endif
686	! endif
687
688	ENDDO
689	ENDDO
690	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
691
692
693	c bilball : uniform albedo, thermal inertia
694	c -----------------------------------------
695	else if (trim(modif) .eq. 'bilball') then
696	write(,) 'constante albedo and iner.therm:'
697	write(,) 'New value for albedo (ex: 0.25) ?'
698	101 read(,,iostat=ierr) alb_bb
699	if(ierr.ne.0) goto 101
700	write(,)
701	write(,) ' uniform albedo (new value):',alb_bb
702	write(,)
703
704	write(,) 'New value for thermal inertia (eg: 247) ?'
705	102 read(,,iostat=ierr) ith_bb
706	if(ierr.ne.0) goto 102
707	write(,) 'uniform thermal inertia (new value):',ith_bb
708	DO j=1,jjp1
709	DO i=1,iip1
710	alb(i,j) = alb_bb ! albedo
711	do isoil=1,nsoilmx
712	ith(i,j,isoil) = ith_bb ! thermal inertia
713	enddo
714	END DO
715	END DO
716	! CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ith,ithfi)
717	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
718	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
719
720	c coldspole : sous-sol de la calotte sud toujours froid
721	c -----------------------------------------------------
722	else if (trim(modif) .eq. 'coldspole') then
723	write(,)'new value for the subsurface temperature',
724	& ' beneath the permanent southern polar cap ? (eg: 141 K)'
725	103 read(,,iostat=ierr) tsud
726	if(ierr.ne.0) goto 103
727	write(,)
728	write(,) ' new value of the subsurface temperature:',tsud
729	c nouvelle temperature sous la calotte permanente
730	do l=2,nsoilmx
731	tsoil(ngridmx,l) = tsud
732	end do
733
734
735	write(,)'new value for the albedo',
736	& 'of the permanent southern polar cap ? (eg: 0.75)'
737	104 read(,,iostat=ierr) albsud
738	if(ierr.ne.0) goto 104
739	write(,)
740
741	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
742	c Option 1: only the albedo of the pole is modified :
743	albfi(ngridmx)=albsud
744	write(,) 'ig=',ngridmx,' albedo perennial cap ',
745	& albfi(ngridmx)
746
747	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
748	c Option 2 A haute resolution : coordonnee de la vrai calotte ~
749	c DO j=1,jjp1
750	c DO i=1,iip1
751	c ig=1+(j-2)*iim +i
752	c if(j.eq.1) ig=1
753	c if(j.eq.jjp1) ig=ngridmx
754	c if ((rlatu(j)*180./pi.lt.-84.).and.
755	c & (rlatu(j)*180./pi.gt.-91.).and.
756	c & (rlonv(i)*180./pi.gt.-91.).and.
757	c & (rlonv(i)*180./pi.lt.0.)) then
758	cc albedo de la calotte permanente fixe a albsud
759	c alb(i,j)=albsud
760	c write(,) 'lat=',rlatu(j)*180./pi,
761	c & ' lon=',rlonv(i)*180./pi
762	cc fin de la condition sur les limites de la calotte permanente
763	c end if
764	c ENDDO
765	c ENDDO
766	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
767
768	c CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
769
770
771	c ptot : Modification of the total pressure: ice + current atmosphere
772	c -------------------------------------------------------------------
773	else if (trim(modif).eq.'ptot') then
774
775	! check if we have a co2_ice surface tracer:
776	if (igcm_co2_ice.eq.0) then
777	write(,) " No surface CO2 ice !"
778	write(,) " only atmospheric pressure will be considered!"
779	endif
780	c calcul de la pression totale glace + atm actuelle
781	patm=0.
782	airetot=0.
783	pcap=0.
784	DO j=1,jjp1
785	DO i=1,iim
786	ig=1+(j-2)*iim +i
787	if(j.eq.1) ig=1
788	if(j.eq.jjp1) ig=ngridmx
789	patm = patm + ps(i,j)*aire(i,j)
790	airetot= airetot + aire(i,j)
791	if (igcm_co2_ice.ne.0) then
792	!pcap = pcap + aire(i,j)co2ice(ig)g
793	pcap = pcap + aire(i,j)qsurf(ig,igcm_co2_ice)g
794	endif
795	ENDDO
796	ENDDO
797	ptoto = pcap + patm
798
799	print*,'Current total pressure at surface (co2 ice + atm) ',
800	& ptoto/airetot
801
802	print*,'new value?'
803	read(,) ptotn
804	ptotn=ptotn*airetot
805	patmn=ptotn-pcap
806	print*,'ptoto,patm,ptotn,patmn'
807	print*,ptoto,patm,ptotn,patmn
808	print*,'Mult. factor for pressure (atm only)', patmn/patm
809	do j=1,jjp1
810	do i=1,iip1
811	ps(i,j)=ps(i,j)*patmn/patm
812	enddo
813	enddo
814
815
816
817	c Correction pour la conservation des traceurs
818	yes=' '
819	do while ((yes.ne.'y').and.(yes.ne.'n'))
820	write(,) 'Do you wish to conserve tracer total mass (y)',
821	& ' or tracer mixing ratio (n) ?'
822	read(*,fmt='(a)') yes
823	end do
824
825	if (yes.eq.'y') then
826	write(,) 'OK : conservation of tracer total mass'
827	DO iq =1, nqtot
828	DO l=1,llm
829	DO j=1,jjp1
830	DO i=1,iip1
831	q(i,j,l,iq)=q(i,j,l,iq)*patm/patmn
832	ENDDO
833	ENDDO
834	ENDDO
835	ENDDO
836	else
837	write(,) 'OK : conservation of tracer mixing ratio'
838	end if
839
840	c Correction pour la pression au niveau de la mer
841	yes=' '
842	do while ((yes.ne.'y').and.(yes.ne.'n'))
843	write(,) 'Do you wish fix pressure at sea level (y)',
844	& ' or not (n) ?'
845	read(*,fmt='(a)') yes
846	end do
847
848	if (yes.eq.'y') then
849	write(,) 'Value?'
850	read(,,iostat=ierr) psea
851	DO i=1,iip1
852	DO j=1,jjp1
853	ps(i,j)=psea
854
855	ENDDO
856	ENDDO
857	write(,) 'psea=',psea
858	else
859	write(,) 'no'
860	end if
861
862
863	c emis : change surface emissivity (added by RW)
864	c ----------------------------------------------
865	else if (trim(modif).eq.'emis') then
866
867	print*,'new value?'
868	read(,) emisread
869
870	do i=1,ngridmx
871	emis(i)=emisread
872	enddo
873
874	c qname : change tracer name
875	c --------------------------
876	else if (trim(modif).eq.'qname') then
877	yes='y'
878	do while (yes.eq.'y')
879	write(,) 'Which tracer name do you want to change ?'
880	do iq=1,nqtot
881	write(*,'(i3,a3,a20)')iq,' : ',trim(tname(iq))
882	enddo
883	write(*,'(a35,i3)')
884	& '(enter tracer number; between 1 and ',nqtot
885	write(,)' or any other value to quit this option)'
886	read(,) iq
887	if ((iq.ge.1).and.(iq.le.nqtot)) then
888	write(,)'Change tracer name ',trim(tname(iq)),' to ?'
889	read(,) txt
890	tname(iq)=txt
891	write(,)'Do you want to change another tracer name (y/n)?'
892	read(*,'(a)') yes
893	else
894	! inapropiate value of iq; quit this option
895	yes='n'
896	endif ! of if ((iq.ge.1).and.(iq.le.nqtot))
897	enddo ! of do while (yes.ne.'y')
898
899	c q=0 : set tracers to zero
900	c -------------------------
901	else if (trim(modif).eq.'q=0') then
902	c mise a 0 des q (traceurs)
903	write(,) 'Tracers set to 0 (1.E-30 in fact)'
904	DO iq =1, nqtot
905	DO l=1,llm
906	DO j=1,jjp1
907	DO i=1,iip1
908	q(i,j,l,iq)=1.e-30
909	ENDDO
910	ENDDO
911	ENDDO
912	ENDDO
913
914	c set surface tracers to zero
915	DO iq =1, nqtot
916	DO ig=1,ngridmx
917	qsurf(ig,iq)=0.
918	ENDDO
919	ENDDO
920
921	c q=x : initialise tracer manually
922	c --------------------------------
923	else if (trim(modif).eq.'q=x') then
924	write(,) 'Which tracer do you want to modify ?'
925	do iq=1,nqtot
926	write(,)iq,' : ',trim(tname(iq))
927	enddo
928	write(,) '(choose between 1 and ',nqtot,')'
929	read(,) iq
930	write(,)'mixing ratio of tracer ',trim(tname(iq)),
931	& ' ? (kg/kg)'
932	read(,) val
933	DO l=1,llm
934	DO j=1,jjp1
935	DO i=1,iip1
936	q(i,j,l,iq)=val
937	ENDDO
938	ENDDO
939	ENDDO
940	write(,) 'SURFACE value of tracer ',trim(tname(iq)),
941	& ' ? (kg/m2)'
942	read(,) val
943	DO ig=1,ngridmx
944	qsurf(ig,iq)=val
945	ENDDO
946
947	c t=profile : initialize temperature with a given profile
948	c -------------------------------------------------------
949	else if (trim(modif) .eq. 't=profile') then
950	write(,) 'Temperature profile from ASCII file'
951	write(,) "'profile.in' e.g. 1D output"
952	write(,) "(one value per line in file; starting with"
953	write(,) "surface value, the 1st atmospheric layer"
954	write(,) "followed by 2nd, etc. up to top of atmosphere)"
955	txt="profile.in"
956	open(33,file=trim(txt),status='old',form='formatted',
957	& iostat=ierr)
958	if (ierr.eq.0) then
959	! OK, found file 'profile_...', load the profile
960	do l=1,llm+1
961	read(33,*,iostat=ierr) profile(l)
962	write(,) profile(l)
963	if (ierr.ne.0) then ! something went wrong
964	exit ! quit loop
965	endif
966	enddo
967	if (ierr.eq.0) then
968	tsurf(1:ngridmx)=profile(1)
969	tsoil(1:ngridmx,1:nsoilmx)=profile(1)
970	do l=1,llm
971	Tset(1:iip1,1:jjp1,l)=profile(l+1)
972	flagtset=.true.
973	enddo
974	ucov(1:iip1,1:jjp1,1:llm)=0.
975	vcov(1:iip1,1:jjm,1:llm)=0.
976	q2(1:ngridmx,1:llm+1)=0.
977	else
978	write(,)'problem reading file ',trim(txt),' !'
979	write(,)'No modifications to temperature'
980	endif
981	else
982	write(,)'Could not find file ',trim(txt),' !'
983	write(,)'No modifications to temperature'
984	endif
985
986	c q=profile : initialize tracer with a given profile
987	c --------------------------------------------------
988	else if (trim(modif) .eq. 'q=profile') then
989	write(,) 'Tracer profile will be sought in ASCII file'
990	write(,) "'profile_tracer' where 'tracer' is tracer name"
991	write(,) "(one value per line in file; starting with"
992	write(,) "surface value, the 1st atmospheric layer"
993	write(,) "followed by 2nd, etc. up to top of atmosphere)"
994	write(,) 'Which tracer do you want to set?'
995	do iq=1,nqtot
996	write(,)iq,' : ',trim(tname(iq))
997	enddo
998	write(,) '(choose between 1 and ',nqtot,')'
999	read(,) iq
1000	if ((iq.lt.1).or.(iq.gt.nqtot)) then
1001	! wrong value for iq, go back to menu
1002	write(,) "wrong input value:",iq
1003	cycle
1004	endif
1005	! look for input file 'profile_tracer'
1006	txt="profile_"//trim(tname(iq))
1007	open(41,file=trim(txt),status='old',form='formatted',
1008	& iostat=ierr)
1009	if (ierr.eq.0) then
1010	! OK, found file 'profile_...', load the profile
1011	do l=1,llm+1
1012	read(41,*,iostat=ierr) profile(l)
1013	if (ierr.ne.0) then ! something went wrong
1014	exit ! quit loop
1015	endif
1016	enddo
1017	if (ierr.eq.0) then
1018	! initialize tracer values
1019	qsurf(:,iq)=profile(1)
1020	do l=1,llm
1021	q(:,:,l,iq)=profile(l+1)
1022	enddo
1023	write(,)'OK, tracer ',trim(tname(iq)),' initialized '
1024	& ,'using values from file ',trim(txt)
1025	else
1026	write(,)'problem reading file ',trim(txt),' !'
1027	write(,)'No modifications to tracer ',trim(tname(iq))
1028	endif
1029	else
1030	write(,)'Could not find file ',trim(txt),' !'
1031	write(,)'No modifications to tracer ',trim(tname(iq))
1032	endif
1033
1034
1035	c wetstart : wet atmosphere with a north to south gradient
1036	c --------------------------------------------------------
1037	else if (trim(modif) .eq. 'wetstart') then
1038	! check that there is indeed a water vapor tracer
1039	if (igcm_h2o_vap.eq.0) then
1040	write(,) "No water vapour tracer! Can't use this option"
1041	stop
1042	endif
1043	DO l=1,llm
1044	DO j=1,jjp1
1045	DO i=1,iip1
1046	q(i,j,l,igcm_h2o_vap)=150.e-6 * (rlatu(j)+pi/2.) / pi
1047	ENDDO
1048	ENDDO
1049	ENDDO
1050
1051	write(,) 'Water mass mixing ratio at north pole='
1052	* ,q(1,1,1,igcm_h2o_vap)
1053	write(,) '---------------------------south pole='
1054	* ,q(1,jjp1,1,igcm_h2o_vap)
1055
1056	c noglacier : remove tropical water ice (to initialize high res sim)
1057	c --------------------------------------------------
1058	else if (trim(modif) .eq. 'noglacier') then
1059	if (igcm_h2o_ice.eq.0) then
1060	write(,) "No water ice tracer! Can't use this option"
1061	stop
1062	endif
1063	do ig=1,ngridmx
1064	j=(ig-2)/iim +2
1065	if(ig.eq.1) j=1
1066	write(,) 'OK: remove surface ice for \|lat\|<45'
1067	if (abs(rlatu(j)*180./pi).lt.45.) then
1068	qsurf(ig,igcm_h2o_ice)=0.
1069	end if
1070	end do
1071
1072
1073	c watercapn : H20 ice on permanent northern cap
1074	c --------------------------------------------------
1075	else if (trim(modif) .eq. 'watercapn') then
1076	if (igcm_h2o_ice.eq.0) then
1077	write(,) "No water ice tracer! Can't use this option"
1078	stop
1079	endif
1080
1081	DO j=1,jjp1
1082	DO i=1,iim
1083	ig=1+(j-2)*iim +i
1084	if(j.eq.1) ig=1
1085	if(j.eq.jjp1) ig=ngridmx
1086
1087	if (rlatu(j)*180./pi.gt.80.) then
1088	qsurf(ig,igcm_h2o_ice)=3.4e3
1089	!do isoil=1,nsoilmx
1090	! ith(i,j,isoil) = 18000. ! thermal inertia
1091	!enddo
1092	write(,)' ==> Ice mesh North boundary (deg)= ',
1093	& rlatv(j-1)*180./pi
1094	end if
1095	ENDDO
1096	ENDDO
1097	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
1098
1099	c$$$ do ig=1,ngridmx
1100	c$$$ j=(ig-2)/iim +2
1101	c$$$ if(ig.eq.1) j=1
1102	c$$$ if (rlatu(j)*180./pi.gt.80.) then
1103	c$$$
1104	c$$$ qsurf(ig,igcm_h2o_ice)=1.e5
1105	c$$$ qsurf(ig,igcm_h2o_vap)=0.0!1.e5
1106	c$$$
1107	c$$$ write(,) 'ig=',ig,' H2O ice mass (kg/m2)= ',
1108	c$$$ & qsurf(ig,igcm_h2o_ice)
1109	c$$$
1110	c$$$ write(,)' ==> Ice mesh South boundary (deg)= ',
1111	c$$$ & rlatv(j)*180./pi
1112	c$$$ end if
1113	c$$$ enddo
1114
1115	c watercaps : H20 ice on permanent southern cap
1116	c -------------------------------------------------
1117	else if (trim(modif) .eq. 'watercaps') then
1118	if (igcm_h2o_ice.eq.0) then
1119	write(,) "No water ice tracer! Can't use this option"
1120	stop
1121	endif
1122
1123	DO j=1,jjp1
1124	DO i=1,iim
1125	ig=1+(j-2)*iim +i
1126	if(j.eq.1) ig=1
1127	if(j.eq.jjp1) ig=ngridmx
1128
1129	if (rlatu(j)*180./pi.lt.-80.) then
1130	qsurf(ig,igcm_h2o_ice)=3.4e3
1131	!do isoil=1,nsoilmx
1132	! ith(i,j,isoil) = 18000. ! thermal inertia
1133	!enddo
1134	write(,)' ==> Ice mesh South boundary (deg)= ',
1135	& rlatv(j-1)*180./pi
1136	end if
1137	ENDDO
1138	ENDDO
1139	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
1140
1141	c$$$ do ig=1,ngridmx
1142	c$$$ j=(ig-2)/iim +2
1143	c$$$ if(ig.eq.1) j=1
1144	c$$$ if (rlatu(j)*180./pi.lt.-80.) then
1145	c$$$ qsurf(ig,igcm_h2o_ice)=1.e5
1146	c$$$ qsurf(ig,igcm_h2o_vap)=0.0 !1.e5
1147	c$$$
1148	c$$$ write(,) 'ig=',ig,' H2O ice mass (kg/m2)= ',
1149	c$$$ & qsurf(ig,igcm_h2o_ice)
1150	c$$$ write(,)' ==> Ice mesh North boundary (deg)= ',
1151	c$$$ & rlatv(j-1)*180./pi
1152	c$$$ end if
1153	c$$$ enddo
1154
1155
1156	c noacglac : H2O ice across highest terrain
1157	c --------------------------------------------
1158	else if (trim(modif) .eq. 'noacglac') then
1159	if (igcm_h2o_ice.eq.0) then
1160	write(,) "No water ice tracer! Can't use this option"
1161	stop
1162	endif
1163	DO j=1,jjp1
1164	DO i=1,iim
1165	ig=1+(j-2)*iim +i
1166	if(j.eq.1) ig=1
1167	if(j.eq.jjp1) ig=ngridmx
1168
1169	if(phis(i,j).gt.1000.*g)then
1170	alb(i,j) = 0.5 ! snow value
1171	do isoil=1,nsoilmx
1172	ith(i,j,isoil) = 18000. ! thermal inertia
1173	! this leads to rnat set to 'ocean' in physiq.F90
1174	! actually no, because it is soil not surface
1175	enddo
1176	endif
1177	ENDDO
1178	ENDDO
1179	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
1180	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
1181	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
1182
1183
1184
1185	c oborealis : H2O oceans across Vastitas Borealis
1186	c -----------------------------------------------
1187	else if (trim(modif) .eq. 'oborealis') then
1188	if (igcm_h2o_ice.eq.0) then
1189	write(,) "No water ice tracer! Can't use this option"
1190	stop
1191	endif
1192	DO j=1,jjp1
1193	DO i=1,iim
1194	ig=1+(j-2)*iim +i
1195	if(j.eq.1) ig=1
1196	if(j.eq.jjp1) ig=ngridmx
1197
1198	if(phis(i,j).lt.-4000.*g)then
1199	! if( (phis(i,j).lt.-4000.*g)
1200	! & .and. (rlatu(j)*180./pi.lt.0.) )then ! south hemisphere only
1201	! phis(i,j)=-8000.0*g ! proper ocean
1202	phis(i,j)=-4000.0*g ! scanty ocean
1203
1204	alb(i,j) = 0.07 ! oceanic value
1205	do isoil=1,nsoilmx
1206	ith(i,j,isoil) = 18000. ! thermal inertia
1207	! this leads to rnat set to 'ocean' in physiq.F90
1208	! actually no, because it is soil not surface
1209	enddo
1210	endif
1211	ENDDO
1212	ENDDO
1213	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
1214	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
1215	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
1216
1217	c iborealis : H2O ice in Northern plains
1218	c --------------------------------------
1219	else if (trim(modif) .eq. 'iborealis') then
1220	if (igcm_h2o_ice.eq.0) then
1221	write(,) "No water ice tracer! Can't use this option"
1222	stop
1223	endif
1224	DO j=1,jjp1
1225	DO i=1,iim
1226	ig=1+(j-2)*iim +i
1227	if(j.eq.1) ig=1
1228	if(j.eq.jjp1) ig=ngridmx
1229
1230	if(phis(i,j).lt.-4000.*g)then
1231	!qsurf(ig,igcm_h2o_ice)=1.e3
1232	qsurf(ig,igcm_h2o_ice)=241.4 ! to make total 33 kg m^-2
1233	endif
1234	ENDDO
1235	ENDDO
1236	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
1237	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
1238	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
1239
1240
1241	c oceanball : H2O liquid everywhere
1242	c ----------------------------
1243	else if (trim(modif) .eq. 'oceanball') then
1244	if (igcm_h2o_ice.eq.0) then
1245	write(,) "No water ice tracer! Can't use this option"
1246	stop
1247	endif
1248	DO j=1,jjp1
1249	DO i=1,iim
1250	ig=1+(j-2)*iim +i
1251	if(j.eq.1) ig=1
1252	if(j.eq.jjp1) ig=ngridmx
1253
1254	qsurf(ig,igcm_h2o_vap)=0.0 ! for ocean, this is infinite source
1255	qsurf(ig,igcm_h2o_ice)=0.0
1256	alb(i,j) = 0.07 ! ocean value
1257
1258	do isoil=1,nsoilmx
1259	ith(i,j,isoil) = 18000. ! thermal inertia
1260	!ith(i,j,isoil) = 50. ! extremely low for test
1261	! this leads to rnat set to 'ocean' in physiq.F90
1262	enddo
1263
1264	ENDDO
1265	ENDDO
1266	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
1267	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
1268	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
1269
1270	c iceball : H2O ice everywhere
1271	c ----------------------------
1272	else if (trim(modif) .eq. 'iceball') then
1273	if (igcm_h2o_ice.eq.0) then
1274	write(,) "No water ice tracer! Can't use this option"
1275	stop
1276	endif
1277	DO j=1,jjp1
1278	DO i=1,iim
1279	ig=1+(j-2)*iim +i
1280	if(j.eq.1) ig=1
1281	if(j.eq.jjp1) ig=ngridmx
1282
1283	qsurf(ig,igcm_h2o_vap)=-50. ! for ocean, this is infinite source
1284	qsurf(ig,igcm_h2o_ice)=50. ! == to 0.05 m of oceanic ice
1285	alb(i,j) = 0.6 ! ice albedo value
1286
1287	do isoil=1,nsoilmx
1288	!ith(i,j,isoil) = 18000. ! thermal inertia
1289	! this leads to rnat set to 'ocean' in physiq.F90
1290	enddo
1291
1292	ENDDO
1293	ENDDO
1294	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
1295	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
1296
1297	c isotherm : Isothermal temperatures and no winds
1298	c -----------------------------------------------
1299	else if (trim(modif) .eq. 'isotherm') then
1300
1301	write(,)'Isothermal temperature of the atmosphere,
1302	& surface and subsurface'
1303	write(,) 'Value of this temperature ? :'
1304	203 read(,,iostat=ierr) Tiso
1305	if(ierr.ne.0) goto 203
1306
1307	tsurf(1:ngridmx)=Tiso
1308
1309	tsoil(1:ngridmx,1:nsoilmx)=Tiso
1310
1311	Tset(1:iip1,1:jjp1,1:llm)=Tiso
1312	flagtset=.true.
1313
1314	ucov(1:iip1,1:jjp1,1:llm)=0
1315	vcov(1:iip1,1:jjm,1:llm)=0
1316	q2(1:ngridmx,1:llm+1)=0
1317
1318	c radequi : Radiative equilibrium profile of temperatures and no winds
1319	c --------------------------------------------------------------------
1320	else if (trim(modif) .eq. 'radequi') then
1321
1322	write(,)'radiative equilibrium temperature profile'
1323
1324	do ig=1, ngridmx
1325	teque= 335.0-60.0sin(latfi(ig))sin(latfi(ig))-
1326	& 10.0cos(latfi(ig))cos(latfi(ig))
1327	tsurf(ig) = MAX(220.0,teque)
1328	end do
1329	do l=2,nsoilmx
1330	do ig=1, ngridmx
1331	tsoil(ig,l) = tsurf(ig)
1332	end do
1333	end do
1334	DO j=1,jjp1
1335	DO i=1,iim
1336	Do l=1,llm
1337	teque=335.-60.0sin(rlatu(j))sin(rlatu(j))-
1338	& 10.0cos(rlatu(j))cos(rlatu(j))
1339	Tset(i,j,l)=MAX(220.0,teque)
1340	end do
1341	end do
1342	end do
1343	flagtset=.true.
1344	ucov(1:iip1,1:jjp1,1:llm)=0
1345	vcov(1:iip1,1:jjm,1:llm)=0
1346	q2(1:ngridmx,1:llm+1)=0
1347
1348	c coldstart : T set 1K above CO2 frost point and no winds
1349	c ------------------------------------------------
1350	else if (trim(modif) .eq. 'coldstart') then
1351
1352	write(,)'set temperature of the atmosphere,'
1353	&,'surface and subsurface how many degrees above CO2 frost point?'
1354	204 read(,,iostat=ierr) Tabove
1355	if(ierr.ne.0) goto 204
1356
1357	DO j=1,jjp1
1358	DO i=1,iim
1359	ig=1+(j-2)*iim +i
1360	if(j.eq.1) ig=1
1361	if(j.eq.jjp1) ig=ngridmx
1362	tsurf(ig) = (-3167.8)/(log(.01*ps(i,j))-23.23)+Tabove
1363	END DO
1364	END DO
1365	do l=1,nsoilmx
1366	do ig=1, ngridmx
1367	tsoil(ig,l) = tsurf(ig)
1368	end do
1369	end do
1370	DO j=1,jjp1
1371	DO i=1,iim
1372	Do l=1,llm
1373	pp = aps(l) +bps(l)*ps(i,j)
1374	Tset(i,j,l)=(-3167.8)/(log(.01*pp)-23.23)+Tabove
1375	end do
1376	end do
1377	end do
1378
1379	flagtset=.true.
1380	ucov(1:iip1,1:jjp1,1:llm)=0
1381	vcov(1:iip1,1:jjm,1:llm)=0
1382	q2(1:ngridmx,1:llm+1)=0
1383
1384
1385	c co2ice=0 : remove CO2 polar ice caps'
1386	c ------------------------------------------------
1387	else if (trim(modif) .eq. 'co2ice=0') then
1388	! check that there is indeed a co2_ice tracer ...
1389	if (igcm_co2_ice.ne.0) then
1390	do ig=1,ngridmx
1391	!co2ice(ig)=0
1392	qsurf(ig,igcm_co2_ice)=0
1393	emis(ig)=emis(ngridmx/2)
1394	end do
1395	else
1396	write(,) "Can't remove CO2 ice!! (no co2_ice tracer)"
1397	endif
1398
1399	! therm_ini_s: (re)-set soil thermal inertia to reference surface values
1400	! ----------------------------------------------------------------------
1401
1402	else if (trim(modif) .eq. 'therm_ini_s') then
1403	! write(,)"surfithfi(1):",surfithfi(1)
1404	do isoil=1,nsoilmx
1405	inertiedat(1:ngridmx,isoil)=surfithfi(1:ngridmx)
1406	enddo
1407	write(,)'OK: Soil thermal inertia has been reset to referenc
1408	&e surface values'
1409	! write(,)"inertiedat(1,1):",inertiedat(1,1)
1410	ithfi(:,:)=inertiedat(:,:)
1411	! recast ithfi() onto ith()
1412	call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
1413	! Check:
1414	! do i=1,iip1
1415	! do j=1,jjp1
1416	! do isoil=1,nsoilmx
1417	! write(77,*) i,j,isoil," ",ith(i,j,isoil)
1418	! enddo
1419	! enddo
1420	! enddo
1421
1422
1423
1424	c slab_ocean_initialisation
1425	c ------------------------------------------------
1426	else if (trim(modif) .eq. 'slab_ocean_0') then
1427	write(,)'OK: initialisation of slab ocean'
1428
1429	DO ig=1, ngridmx
1430	rnat(ig)=1.
1431	tslab(ig,1)=0.
1432	tslab(ig,2)=0.
1433	tsea_ice(ig)=0.
1434	sea_ice(ig)=0.
1435	pctsrf_sic(ig)=0.
1436
1437	if(ithfi(ig,1).GT.10000.)then
1438	rnat(ig)=0.
1439	phisfi(ig)=0.
1440	tsea_ice(ig)=273.15-1.8
1441	tslab(ig,1)=tsurf(ig)
1442	tslab(ig,2)=tsurf(ig)!*2./3.+(273.15-1.8)/3.
1443	endif
1444
1445	ENDDO
1446	CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
1447
1448
1449
1450	else
1451	write(,) ' Unknown (misspelled?) option!!!'
1452	end if ! of if (trim(modif) .eq. '...') elseif ...
1453
1454
1455
1456	enddo ! of do ! infinite loop on liste of changes
1457
1458	999 continue
1459
1460
1461	c=======================================================================
1462	c Initialisation for cloud fraction and oceanic ice (added by BC 2010)
1463	c=======================================================================
1464	DO ig=1, ngridmx
1465	totalfrac(ig)=0.5
1466	DO l=1,llm
1467	cloudfrac(ig,l)=0.5
1468	ENDDO
1469	! Ehouarn, march 2012: also add some initialisation for hice
1470	hice(ig)=0.0
1471	ENDDO
1472
1473	c========================================================
1474
1475	! DO ig=1,ngridmx
1476	! IF(tsurf(ig) .LT. 273.16-1.8) THEN
1477	! hice(ig)=(273.16-1.8-tsurf(ig))/(273.16-1.8-240)*1.
1478	! hice(ig)=min(hice(ig),1.0)
1479	! ENDIF
1480	! ENDDO
1481
1482
1483
1484
1485	c=======================================================================
1486	c Correct pressure on the new grid (menu 0)
1487	c=======================================================================
1488
1489
1490	if ((choix_1.eq.0).and.(.not.flagps0)) then
1491	r = 1000.*8.31/mugaz
1492
1493	do j=1,jjp1
1494	do i=1,iip1
1495	ps(i,j) = ps(i,j) *
1496	. exp((phisold_newgrid(i,j)-phis(i,j)) /
1497	. (t(i,j,1) * r))
1498	end do
1499	end do
1500
1501	c periodicite de ps en longitude
1502	do j=1,jjp1
1503	ps(1,j) = ps(iip1,j)
1504	end do
1505	end if
1506
1507
1508	c=======================================================================
1509	c=======================================================================
1510
1511	c=======================================================================
1512	c Initialisation de la physique / ecriture de newstartfi :
1513	c=======================================================================
1514
1515
1516	CALL inifilr
1517	CALL pression(ip1jmp1, ap, bp, ps, p3d)
1518
1519	c-----------------------------------------------------------------------
1520	c Initialisation pks:
1521	c-----------------------------------------------------------------------
1522
1523	CALL exner_hyb(ip1jmp1, ps, p3d, beta, pks, pk, pkf)
1524	! Calcul de la temperature potentielle teta
1525
1526	if (flagtset) then
1527	DO l=1,llm
1528	DO j=1,jjp1
1529	DO i=1,iim
1530	teta(i,j,l) = Tset(i,j,l) * cpp/pk(i,j,l)
1531	ENDDO
1532	teta (iip1,j,l)= teta (1,j,l)
1533	ENDDO
1534	ENDDO
1535	else if (choix_1.eq.0) then
1536	DO l=1,llm
1537	DO j=1,jjp1
1538	DO i=1,iim
1539	teta(i,j,l) = t(i,j,l) * cpp/pk(i,j,l)
1540	ENDDO
1541	teta (iip1,j,l)= teta (1,j,l)
1542	ENDDO
1543	ENDDO
1544	endif
1545
1546	C Calcul intermediaire
1547	c
1548	if (choix_1.eq.0) then
1549	CALL massdair( p3d, masse )
1550	c
1551	print *,' ALPHAX ',alphax
1552
1553	DO l = 1, llm
1554	DO i = 1, iim
1555	xppn(i) = aire( i, 1 ) * masse( i , 1 , l )
1556	xpps(i) = aire( i,jjp1 ) * masse( i , jjp1 , l )
1557	ENDDO
1558	xpn = SUM(xppn)/apoln
1559	xps = SUM(xpps)/apols
1560	DO i = 1, iip1
1561	masse( i , 1 , l ) = xpn
1562	masse( i , jjp1 , l ) = xps
1563	ENDDO
1564	ENDDO
1565	endif
1566	phis(iip1,:) = phis(1,:)
1567
1568	itau=0
1569	if (choix_1.eq.0) then
1570	day_ini=int(date)
1571	endif
1572	c
1573	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
1574
1575	CALL caldyn0( itau,ucov,vcov,teta,ps,masse,pk,phis ,
1576	* phi,w, pbaru,pbarv,day_ini+time )
1577
1578
1579	CALL dynredem0("restart.nc",day_ini,anneeref,phis,nqtot)
1580	CALL dynredem1("restart.nc",0.0,vcov,ucov,teta,q,nqtot,masse,ps)
1581	C
1582	C Ecriture etat initial physique
1583	C
1584
1585	call physdem0("restartfi.nc",lonfi,latfi,nsoilmx,ngridmx,llm,
1586	& nqtot,dtphys,real(day_ini),0.0,
1587	& airefi,albfi,ithfi,zmea,zstd,zsig,zgam,zthe)
1588	call physdem1("restartfi.nc",nsoilmx,ngridmx,llm,nq,
1589	& dtphys,real(day_ini),
1590	& tsurf,tsoil,emis,q2,qsurf,
1591	& cloudfrac,totalfrac,hice,
1592	& rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
1593
1594	c=======================================================================
1595	c Formats
1596	c=======================================================================
1597
1598	1 FORMAT(//10x,'la valeur de im =',i4,2x,'lue sur le fichier de dema
1599	*rrage est differente de la valeur parametree iim =',i4//)
1600	2 FORMAT(//10x,'la valeur de jm =',i4,2x,'lue sur le fichier de dema
1601	*rrage est differente de la valeur parametree jjm =',i4//)
1602	3 FORMAT(//10x,'la valeur de lllm =',i4,2x,'lue sur le fichier demar
1603	*rage est differente de la valeur parametree llm =',i4//)
1604
1605	write(,) "newstart: All is well that ends well."
1606
1607	end
1608

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

Download in other formats:

Original Format