Context Navigation

newstart.F @ 1303

Last change on this file since 1303 was 1297, checked in by bclmd, 10 years ago
LMDZ.GENERIC: slab ocean added
File size: 54.7 KB

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

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

Download in other formats:

Original Format