Context Navigation

newstart.F @ 937

Last change on this file since 937 was 905, checked in by emillour, 12 years ago

Generic GCM:

removed the "-static" ifort compilation option (problematic on Gnome cluster) and added some additional debug options
fixed bug in newstart about initialization of albedo and thermal inertia.

File size: 49.5 KB

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

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

Download in other formats:

Original Format