Context Navigation

newstart.F @ 401

Last change on this file since 401 was 374, checked in by emillour, 14 years ago
Generic GCM: Upgrade: The location of the 'datagcm' directory can now be given in the callphys.def file ( datadir = /absolute/path/to/datagcm ). Changed "datafile.h" into a F90 module "datafile_mod.F90" and spread this change to all routines that used to use "datafile.h". EM
File size: 65.6 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	implicit none
18
19	#include "dimensions.h"
20	#include "dimphys.h"
21	#include "surfdat.h"
22	#include "comsoil.h"
23	#include "planete.h"
24	#include "paramet.h"
25	#include "comconst.h"
26	#include "comvert.h"
27	#include "comgeom2.h"
28	#include "control.h"
29	#include "logic.h"
30	#include "description.h"
31	#include "ener.h"
32	#include "temps.h"
33	#include "lmdstd.h"
34	#include "comdissnew.h"
35	#include "clesph0.h"
36	#include "serre.h"
37	#include "netcdf.inc"
38	#include "advtrac.h"
39	#include "tracer.h"
40	c=======================================================================
41	c Declarations
42	c=======================================================================
43
44	c Variables dimension du fichier "start_archive"
45	c------------------------------------
46	CHARACTER relief*3
47
48
49	c Variables pour les lectures NetCDF des fichiers "start_archive"
50	c--------------------------------------------------
51	INTEGER nid_dyn, nid_fi,nid,nvarid
52	INTEGER length
53	parameter (length = 100)
54	INTEGER tab0
55	INTEGER NB_ETATMAX
56	parameter (NB_ETATMAX = 100)
57
58	REAL date
59	REAL p_rad,p_omeg,p_g,p_cpp,p_mugaz,p_daysec
60
61	c Variable histoire
62	c------------------
63	REAL vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm) ! vents covariants
64	REAL phis(iip1,jjp1)
65	REAL q(iip1,jjp1,llm,nqmx) ! champs advectes
66
67	c autre variables dynamique nouvelle grille
68	c------------------------------------------
69	REAL pks(iip1,jjp1)
70	REAL w(iip1,jjp1,llm+1)
71	REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
72	! REAL dv(ip1jm,llm),du(ip1jmp1,llm)
73	! REAL dh(ip1jmp1,llm),dp(ip1jmp1)
74	REAL phi(iip1,jjp1,llm)
75
76	integer klatdat,klongdat
77	PARAMETER (klatdat=180,klongdat=360)
78
79	c Physique sur grille scalaire
80	c----------------------------
81	real zmeaS(iip1,jjp1),zstdS(iip1,jjp1)
82	real zsigS(iip1,jjp1),zgamS(iip1,jjp1),ztheS(iip1,jjp1)
83
84	c variable physique
85	c------------------
86	REAL tsurf(ngridmx) ! surface temperature
87	REAL tsoil(ngridmx,nsoilmx) ! soil temperature
88	! REAL co2ice(ngridmx) ! CO2 ice layer
89	REAL emis(ngridmx) ! surface emissivity
90	real emisread ! added by RW
91	REAL qsurf(ngridmx,nqmx)
92	REAL q2(ngridmx,nlayermx+1)
93	! REAL rnaturfi(ngridmx)
94	real alb(iip1,jjp1),albfi(ngridmx) ! albedos
95	real ith(iip1,jjp1,nsoilmx),ithfi(ngridmx,nsoilmx) ! thermal inertia (3D)
96	real surfith(iip1,jjp1),surfithfi(ngridmx) ! surface thermal inertia (2D)
97	REAL latfi(ngridmx),lonfi(ngridmx),airefi(ngridmx)
98
99	INTEGER i,j,l,isoil,ig,idum
100	real mugaz ! molar mass of the atmosphere
101
102	integer ierr
103
104	c Variables on the new grid along scalar points
105	c------------------------------------------------------
106	! REAL p(iip1,jjp1)
107	REAL t(iip1,jjp1,llm)
108	REAL tset(iip1,jjp1,llm)
109	real phisold_newgrid(iip1,jjp1)
110	REAL :: teta(iip1, jjp1, llm)
111	REAL :: pk(iip1,jjp1,llm)
112	REAL :: pkf(iip1,jjp1,llm)
113	REAL :: ps(iip1, jjp1)
114	REAL :: masse(iip1,jjp1,llm)
115	REAL :: xpn,xps,xppn(iim),xpps(iim)
116	REAL :: p3d(iip1, jjp1, llm+1)
117	REAL :: beta(iip1,jjp1,llm)
118	! REAL dteta(ip1jmp1,llm)
119
120	c Variable de l'ancienne grille
121	c------------------------------
122	real time
123	real tab_cntrl(100)
124	real tab_cntrl_bis(100)
125
126	c variables diverses
127	c-------------------
128	real choix_1,pp
129	character*80 fichnom
130	integer Lmodif,iq,thermo
131	character modif*20
132	real z_reel(iip1,jjp1)
133	real tsud,albsud,alb_bb,ith_bb,Tiso,Tabove
134	real ptoto,pcap,patm,airetot,ptotn,patmn,psea
135	! real ssum
136	character*1 yes
137	logical :: flagtset=.false. , flagps0=.false.
138	real val, val2, val3 ! to store temporary variables
139	real :: iceith=2000 ! thermal inertia of subterranean ice
140	integer iref,jref
141
142	INTEGER :: itau
143
144	INTEGER :: nq,numvanle
145	character(len=20) :: txt ! to store some text
146	integer :: count
147
148	! MONS data:
149	real :: MONS_Hdn(iip1,jjp1) ! Hdn: %WEH=Mass fraction of H2O
150	real :: MONS_d21(iip1,jjp1) ! ice table "depth" (in kg/m2)
151	! coefficient to apply to convert d21 to 'true' depth (m)
152	real :: MONS_coeff
153	real :: MONS_coeffS ! coeff for southern hemisphere
154	real :: MONS_coeffN ! coeff for northern hemisphere
155	! real,parameter :: icedepthmin=1.e-3 ! Ice begins at most at that depth
156
157	! added by RW for test
158	real pmean, phi0
159
160	! added by BC for equilibrium temperature startup
161	real teque
162
163	! added by BC for cloud fraction setup
164	REAL hice(ngridmx),cloudfrac(ngridmx,nlayermx)
165	REAL totalfrac(ngridmx)
166
167	! added by RW for nuketharsis
168	real fact1
169	real fact2
170
171
172	c sortie visu pour les champs dynamiques
173	c---------------------------------------
174	! INTEGER :: visuid
175	! real :: time_step,t_ops,t_wrt
176	! CHARACTER*80 :: visu_file
177
178	cpp = 744.499 ! for Mars, instead of 1004.70885 (Earth)
179	preff = 610. ! for Mars, instead of 101325. (Earth)
180	pa = 20 ! for Mars, instead of 500 (Earth)
181
182	c=======================================================================
183	c Choice of the start file(s) to use
184	c=======================================================================
185
186	write(,) 'From which kind of files do you want to create new',
187	. 'start and startfi files'
188	write(,) ' 0 - from a file start_archive'
189	write(,) ' 1 - from files start and startfi'
190
191	c-----------------------------------------------------------------------
192	c Open file(s) to modify (start or start_archive)
193	c-----------------------------------------------------------------------
194
195	DO
196	read(,,iostat=ierr) choix_1
197	if ((choix_1 /= 0).OR.(choix_1 /=1)) EXIT
198	ENDDO
199
200	c Open start_archive
201	c ~~~~~~~~~~~~~~~~~~~~~~~~~~
202	if (choix_1.eq.0) then
203
204	write(,) 'Creating start files from:'
205	write(,) './start_archive.nc'
206	write(,)
207	fichnom = 'start_archive.nc'
208	ierr = NF_OPEN (fichnom, NF_NOWRITE,nid)
209	IF (ierr.NE.NF_NOERR) THEN
210	write(6,*)' Problem opening file:',fichnom
211	write(6,*)' ierr = ', ierr
212	CALL ABORT
213	ENDIF
214	tab0 = 50
215	Lmodif = 1
216
217	c OR open start and startfi files
218	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
219	else
220	write(,) 'Creating start files from:'
221	write(,) './start.nc and ./startfi.nc'
222	write(,)
223	fichnom = 'start.nc'
224	ierr = NF_OPEN (fichnom, NF_NOWRITE,nid_dyn)
225	IF (ierr.NE.NF_NOERR) THEN
226	write(6,*)' Problem opening file:',fichnom
227	write(6,*)' ierr = ', ierr
228	CALL ABORT
229	ENDIF
230
231	fichnom = 'startfi.nc'
232	ierr = NF_OPEN (fichnom, NF_NOWRITE,nid_fi)
233	IF (ierr.NE.NF_NOERR) THEN
234	write(6,*)' Problem opening file:',fichnom
235	write(6,*)' ierr = ', ierr
236	CALL ABORT
237	ENDIF
238
239	tab0 = 0
240	Lmodif = 0
241
242	endif
243
244	c-----------------------------------------------------------------------
245	c Lecture du tableau des parametres du run (pour la dynamique)
246	c-----------------------------------------------------------------------
247
248	if (choix_1.eq.0) then
249
250	write(,) 'reading tab_cntrl START_ARCHIVE'
251	c
252	ierr = NF_INQ_VARID (nid, "controle", nvarid)
253	#ifdef NC_DOUBLE
254	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tab_cntrl)
255	#else
256	ierr = NF_GET_VAR_REAL(nid, nvarid, tab_cntrl)
257	#endif
258	c
259	else if (choix_1.eq.1) then
260
261	write(,) 'reading tab_cntrl START'
262	c
263	ierr = NF_INQ_VARID (nid_dyn, "controle", nvarid)
264	#ifdef NC_DOUBLE
265	ierr = NF_GET_VAR_DOUBLE(nid_dyn, nvarid, tab_cntrl)
266	#else
267	ierr = NF_GET_VAR_REAL(nid_dyn, nvarid, tab_cntrl)
268	#endif
269	c
270	write(,) 'reading tab_cntrl STARTFI'
271	c
272	ierr = NF_INQ_VARID (nid_fi, "controle", nvarid)
273	#ifdef NC_DOUBLE
274	ierr = NF_GET_VAR_DOUBLE(nid_fi, nvarid, tab_cntrl_bis)
275	#else
276	ierr = NF_GET_VAR_REAL(nid_fi, nvarid, tab_cntrl_bis)
277	#endif
278	c
279	do i=1,50
280	tab_cntrl(i+50)=tab_cntrl_bis(i)
281	enddo
282	write(,) 'printing tab_cntrl', tab_cntrl
283	do i=1,100
284	write(,) i,tab_cntrl(i)
285	enddo
286
287	endif
288	c-----------------------------------------------------------------------
289	c Initialisation des constantes dynamique
290	c-----------------------------------------------------------------------
291
292	kappa = tab_cntrl(9)
293	etot0 = tab_cntrl(12)
294	ptot0 = tab_cntrl(13)
295	ztot0 = tab_cntrl(14)
296	stot0 = tab_cntrl(15)
297	ang0 = tab_cntrl(16)
298	write(,) "Newstart: kappa,etot0,ptot0,ztot0,stot0,ang0"
299	write(,) kappa,etot0,ptot0,ztot0,stot0,ang0
300
301	! for vertical coordinate
302	preff=tab_cntrl(18) ! reference surface pressure
303	pa=tab_cntrl(17) ! reference pressure at which coord is purely pressure
304	!NB: in start_archive files tab_cntrl(17)=tab_cntrl(18)=0
305	if (preff.eq.0) then
306	preff=610
307	pa=20
308	endif
309	write(,) "Newstart: preff=",preff," pa=",pa
310	yes=' '
311	do while ((yes.ne.'y').and.(yes.ne.'n'))
312	write(,) "Change the values of preff and pa? (y/n)"
313	read(*,fmt='(a)') yes
314	end do
315	if (yes.eq.'y') then
316	write(,)"New value of reference surface pressure preff?"
317	write(,)" (for Mars, typically preff=610)"
318	read(,) preff
319	write(,)"New value of reference pressure pa for purely"
320	write(,)"pressure levels (for hybrid coordinates)?"
321	write(,)" (for Mars, typically pa=20)"
322	read(,) pa
323	endif
324	c-----------------------------------------------------------------------
325	c Lecture du tab_cntrl et initialisation des constantes physiques
326	c - pour start: Lmodif = 0 => pas de modifications possibles
327	c (modif dans le tabfi de readfi + loin)
328	c - pour start_archive: Lmodif = 1 => modifications possibles
329	c-----------------------------------------------------------------------
330	if (choix_1.eq.0) then
331	call tabfi (nid,Lmodif,tab0,day_ini,lllm,p_rad,
332	. p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
333	else if (choix_1.eq.1) then
334	call tabfi (nid_fi,Lmodif,tab0,day_ini,lllm,p_rad,
335	. p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
336	endif
337
338	rad = p_rad
339	omeg = p_omeg
340	g = p_g
341	cpp = p_cpp
342	mugaz = p_mugaz
343	daysec = p_daysec
344
345	! print*,'daysec=',p_daysec
346	! stop
347
348	! write(,) 'aire',aire
349
350
351	c=======================================================================
352	c INITIALISATIONS DIVERSES
353	c=======================================================================
354	! Load tracer names:
355	call iniadvtrac(nq,numvanle)
356	! tnom(:) now contains tracer names
357	! Initialize global tracer indexes (stored in tracer.h)
358	call initracer()
359	! Load parameters from run.def file
360	CALL defrun_new( 99, .TRUE. )
361	CALL iniconst
362	CALL inigeom
363	idum=-1
364	idum=0
365
366	c Initialisation coordonnees /aires
367	c -------------------------------
368	! Note: rlatu(:) and rlonv(:) are commons defined in "comgeom.h"
369	! rlatu() and rlonv() are given in radians
370	latfi(1)=rlatu(1)
371	lonfi(1)=0.
372	DO j=2,jjm
373	DO i=1,iim
374	latfi((j-2)*iim+1+i)=rlatu(j)
375	lonfi((j-2)*iim+1+i)=rlonv(i)
376	ENDDO
377	ENDDO
378	latfi(ngridmx)=rlatu(jjp1)
379	lonfi(ngridmx)=0.
380	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
381
382	c=======================================================================
383	c lecture topographie, albedo, inertie thermique, relief sous-maille
384	c=======================================================================
385
386	if (choix_1.ne.1) then ! pour ne pas avoir besoin du fichier
387	! surface.dat dans le cas des start
388
389	c do while((relief(1:3).ne.'mol').AND.(relief(1:3).ne.'pla'))
390	c write(,)
391	c write(,) 'choix du relief (mola,pla)'
392	c write(,) '(Topographie MGS MOLA, plat)'
393	c read(*,fmt='(a3)') relief
394	relief="mola"
395	c enddo
396
397	CALL datareadnc(relief,phis,alb,surfith,zmeaS,zstdS,zsigS,zgamS,
398	. ztheS)
399
400	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
401	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,surfith,surfithfi)
402	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
403	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zmeaS,zmea)
404	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zstdS,zstd)
405	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zsigS,zsig)
406	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zgamS,zgam)
407	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ztheS,zthe)
408
409	endif ! of if (choix_1.ne.1)
410
411
412	c=======================================================================
413	c Lecture des fichiers (start ou start_archive)
414	c=======================================================================
415
416	if (choix_1.eq.0) then
417
418	write(,) 'Reading file START_ARCHIVE'
419	CALL lect_start_archive(date,tsurf,tsoil,emis,q2,
420	. t,ucov,vcov,ps,teta,phisold_newgrid,q,qsurf,
421	& surfith,nid)
422	write(,) "OK, read start_archive file"
423	! copy soil thermal inertia
424	ithfi(:,:)=inertiedat(:,:)
425
426	ierr= NF_CLOSE(nid)
427
428	else if (choix_1.eq.1) then ! c'est l'appel a tabfi de phyeta0 qui
429	! permet de changer les valeurs du
430	! tab_cntrl Lmodif=1
431	tab0=0
432	Lmodif=1 ! Lmodif set to 1 to allow modifications in phyeta0
433	write(,) 'Reading file START'
434	fichnom = 'start.nc'
435	CALL dynetat0(fichnom,nqmx,vcov,ucov,teta,q,masse,
436	. ps,phis,time)
437
438	write(,) 'Reading file STARTFI'
439	fichnom = 'startfi.nc'
440	CALL phyetat0 (fichnom,tab0,Lmodif,nsoilmx,nqmx,
441	. day_ini,time,
442	. tsurf,tsoil,emis,q2,qsurf, !) ! temporary modif by RDW
443	. cloudfrac,totalfrac,hice)
444
445	! copy albedo and soil thermal inertia
446	do i=1,ngridmx
447	albfi(i) = albedodat(i)
448	do j=1,nsoilmx
449	ithfi(i,j) = inertiedat(i,j)
450	enddo
451	! build a surfithfi(:) using 1st layer of ithfi(:), which might
452	! be neede later on if reinitializing soil thermal inertia
453	surfithfi(i)=ithfi(i,1)
454	enddo
455
456	else
457	CALL exit(1)
458	endif
459
460	dtvr = daysec/FLOAT(day_step)
461	dtphys = dtvr * FLOAT(iphysiq)
462
463	c=======================================================================
464	c
465	c=======================================================================
466
467	do ! infinite loop on list of changes
468
469	write(,)
470	write(,)
471	write(,) 'List of possible changes :'
472	write(,) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
473	write(,)
474	write(,) 'flat : no topography ("aquaplanet")'
475	write(,) 'nuketharsis : no Tharsis bulge'
476	write(,) 'bilball : uniform albedo and thermal inertia'
477	write(,) 'coldspole : cold subsurface and high albedo at S.pole'
478	write(,) 'qname : change tracer name'
479	write(,) 'q=0 : ALL tracer =zero'
480	write(,) 'q=x : give a specific uniform value to one tracer'
481	write(,) 'ini_q : tracers initialisation for chemistry, water an
482	$d ice '
483	write(,) 'ini_q-H2O : tracers initialisation for chemistry and
484	$ice '
485	write(,) 'ini_q-iceH2O : tracers initialisation for chemistry on
486	$ly '
487	write(,) 'noglacier : Remove tropical H2O ice if \|lat\|<45'
488	write(,) 'watercapn : H20 ice on permanent N polar cap '
489	write(,) 'watercaps : H20 ice on permanent S polar cap '
490	write(,) 'noacglac : H2O ice across Noachis Terra'
491	write(,) 'oborealis : H2O ice across Vastitas Borealis'
492	write(,) 'iceball : Thick ice layer all over surface'
493	write(,) 'wetstart : start with a wet atmosphere'
494	write(,) 'isotherm : Isothermal Temperatures, wind set to zero'
495	write(,) 'radequi : Earth-like radiative equilibrium temperature
496	$ profile (lat-alt) and winds set to zero'
497	write(,) 'coldstart : Start X K above the CO2 frost point and
498	$set wind to zero (assumes 100% CO2)'
499	write(,) 'co2ice=0 : remove CO2 polar cap'
500	write(,) 'ptot : change total pressure'
501	write(,) 'emis : change surface emissivity'
502	write(,) 'therm_ini_s : Set soil thermal inertia to reference sur
503	&face values'
504	! write(,) 'subsoilice_n : Put deep underground ice layer in northe
505	! &rn hemisphere'
506	! write(,) 'subsoilice_s : Put deep underground ice layer in southe
507	! &rn hemisphere'
508	! write(,) 'mons_ice : Put underground ice layer according to MONS-
509	! &derived data'
510
511	write(,)
512	write(,) 'Change to perform ?'
513	write(,) ' (enter keyword or return to end)'
514	write(,)
515
516	read(*,fmt='(a20)') modif
517	if (modif(1:1) .eq. ' ') exit ! exit loop on changes
518
519	write(,)
520	write(,) trim(modif) , ' : '
521
522	c 'flat : no topography ("aquaplanet")'
523	c -------------------------------------
524	if (modif(1:len_trim(modif)) .eq. 'flat') then
525	c set topo to zero
526	CALL initial0(ip1jmp1,z_reel)
527	CALL multscal(ip1jmp1,z_reel,g,phis)
528	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
529	write(,) 'topography set to zero.'
530	write(,) 'WARNING : the subgrid topography parameters',
531	& ' were not set to zero ! => set calllott to F'
532
533	c Choice of surface pressure
534	yes=' '
535	do while ((yes.ne.'y').and.(yes.ne.'n'))
536	write(,) 'Do you wish to choose homogeneous surface',
537	& 'pressure (y) or let newstart interpolate ',
538	& ' the previous field (n)?'
539	read(*,fmt='(a)') yes
540	end do
541	if (yes.eq.'y') then
542	flagps0=.true.
543	write(,) 'New value for ps (Pa) ?'
544	201 read(,,iostat=ierr) patm
545	if(ierr.ne.0) goto 201
546	write(,)
547	write(,) ' new ps everywhere (Pa) = ', patm
548	write(,)
549	do j=1,jjp1
550	do i=1,iip1
551	ps(i,j)=patm
552	enddo
553	enddo
554	end if
555
556	c 'nuketharsis : no tharsis bulge for Early Mars'
557	c ---------------------------------------------
558	else if (modif(1:len_trim(modif)) .eq. 'nuketharsis') then
559
560	DO j=1,jjp1
561	DO i=1,iim
562	ig=1+(j-2)*iim +i
563	if(j.eq.1) ig=1
564	if(j.eq.jjp1) ig=ngridmx
565
566	fact1=(((rlonv(i)180./pi)+100)*2 +
567	& (rlatu(j)180./pi)2)/65*2
568	fact2=exp( -fact1**2.5 )
569
570	phis(i,j) = phis(i,j) - (phis(i,j)+4000.g)fact2
571
572	! if(phis(i,j).gt.2500.*g)then
573	! if(rlatu(j)*180./pi.gt.-80.)then ! avoid chopping south polar cap
574	! phis(i,j)=2500.*g
575	! endif
576	! endif
577
578	ENDDO
579	ENDDO
580	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
581
582
583	c bilball : uniform albedo, thermal inertia
584	c -----------------------------------------
585	else if (modif(1:len_trim(modif)) .eq. 'bilball') then
586	write(,) 'constante albedo and iner.therm:'
587	write(,) 'New value for albedo (ex: 0.25) ?'
588	101 read(,,iostat=ierr) alb_bb
589	if(ierr.ne.0) goto 101
590	write(,)
591	write(,) ' uniform albedo (new value):',alb_bb
592	write(,)
593
594	write(,) 'New value for thermal inertia (eg: 247) ?'
595	102 read(,,iostat=ierr) ith_bb
596	if(ierr.ne.0) goto 102
597	write(,) 'uniform thermal inertia (new value):',ith_bb
598	DO j=1,jjp1
599	DO i=1,iip1
600	alb(i,j) = alb_bb ! albedo
601	do isoil=1,nsoilmx
602	ith(i,j,isoil) = ith_bb ! thermal inertia
603	enddo
604	END DO
605	END DO
606	! CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ith,ithfi)
607	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
608	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
609
610	c coldspole : sous-sol de la calotte sud toujours froid
611	c -----------------------------------------------------
612	else if (modif(1:len_trim(modif)) .eq. 'coldspole') then
613	write(,)'new value for the subsurface temperature',
614	& ' beneath the permanent southern polar cap ? (eg: 141 K)'
615	103 read(,,iostat=ierr) tsud
616	if(ierr.ne.0) goto 103
617	write(,)
618	write(,) ' new value of the subsurface temperature:',tsud
619	c nouvelle temperature sous la calotte permanente
620	do l=2,nsoilmx
621	tsoil(ngridmx,l) = tsud
622	end do
623
624
625	write(,)'new value for the albedo',
626	& 'of the permanent southern polar cap ? (eg: 0.75)'
627	104 read(,,iostat=ierr) albsud
628	if(ierr.ne.0) goto 104
629	write(,)
630
631	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
632	c Option 1: only the albedo of the pole is modified :
633	albfi(ngridmx)=albsud
634	write(,) 'ig=',ngridmx,' albedo perennial cap ',
635	& albfi(ngridmx)
636
637	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
638	c Option 2 A haute resolution : coordonnee de la vrai calotte ~
639	c DO j=1,jjp1
640	c DO i=1,iip1
641	c ig=1+(j-2)*iim +i
642	c if(j.eq.1) ig=1
643	c if(j.eq.jjp1) ig=ngridmx
644	c if ((rlatu(j)*180./pi.lt.-84.).and.
645	c & (rlatu(j)*180./pi.gt.-91.).and.
646	c & (rlonv(i)*180./pi.gt.-91.).and.
647	c & (rlonv(i)*180./pi.lt.0.)) then
648	cc albedo de la calotte permanente fixe a albsud
649	c alb(i,j)=albsud
650	c write(,) 'lat=',rlatu(j)*180./pi,
651	c & ' lon=',rlonv(i)*180./pi
652	cc fin de la condition sur les limites de la calotte permanente
653	c end if
654	c ENDDO
655	c ENDDO
656	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
657
658	c CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
659
660
661	c ptot : Modification of the total pressure: ice + current atmosphere
662	c -------------------------------------------------------------------
663	else if (modif(1:len_trim(modif)).eq.'ptot') then
664
665	! check if we have a co2_ice surface tracer:
666	if (igcm_co2_ice.eq.0) then
667	write(,) " No surface CO2 ice !"
668	write(,) " only atmospheric pressure will be considered!"
669	endif
670	c calcul de la pression totale glace + atm actuelle
671	patm=0.
672	airetot=0.
673	pcap=0.
674	DO j=1,jjp1
675	DO i=1,iim
676	ig=1+(j-2)*iim +i
677	if(j.eq.1) ig=1
678	if(j.eq.jjp1) ig=ngridmx
679	patm = patm + ps(i,j)*aire(i,j)
680	airetot= airetot + aire(i,j)
681	if (igcm_co2_ice.ne.0) then
682	!pcap = pcap + aire(i,j)co2ice(ig)g
683	pcap = pcap + aire(i,j)qsurf(ig,igcm_co2_ice)g
684	endif
685	ENDDO
686	ENDDO
687	ptoto = pcap + patm
688
689	print*,'Current total pressure at surface (co2 ice + atm) ',
690	& ptoto/airetot
691
692	print*,'new value?'
693	read(,) ptotn
694	ptotn=ptotn*airetot
695	patmn=ptotn-pcap
696	print*,'ptoto,patm,ptotn,patmn'
697	print*,ptoto,patm,ptotn,patmn
698	print*,'Mult. factor for pressure (atm only)', patmn/patm
699	do j=1,jjp1
700	do i=1,iip1
701	ps(i,j)=ps(i,j)*patmn/patm
702	enddo
703	enddo
704
705
706
707	c Correction pour la conservation des traceurs
708	yes=' '
709	do while ((yes.ne.'y').and.(yes.ne.'n'))
710	write(,) 'Do you wish to conserve tracer total mass (y)',
711	& ' or tracer mixing ratio (n) ?'
712	read(*,fmt='(a)') yes
713	end do
714
715	if (yes.eq.'y') then
716	write(,) 'OK : conservation of tracer total mass'
717	DO iq =1, nqmx
718	DO l=1,llm
719	DO j=1,jjp1
720	DO i=1,iip1
721	q(i,j,l,iq)=q(i,j,l,iq)*patm/patmn
722	ENDDO
723	ENDDO
724	ENDDO
725	ENDDO
726	else
727	write(,) 'OK : conservation of tracer mixing ratio'
728	end if
729
730	c Correction pour la pression au niveau de la mer
731	yes=' '
732	do while ((yes.ne.'y').and.(yes.ne.'n'))
733	write(,) 'Do you wish fix pressure at sea level (y)',
734	& ' or not (n) ?'
735	read(*,fmt='(a)') yes
736	end do
737
738	if (yes.eq.'y') then
739	write(,) 'Value?'
740	read(,,iostat=ierr) psea
741	DO i=1,iip1
742	DO j=1,jjp1
743	ps(i,j)=psea
744
745	ENDDO
746	ENDDO
747	write(,) 'psea=',psea
748	else
749	write(,) 'no'
750	end if
751
752
753	c emis : change surface emissivity (added by RW)
754	c ----------------------------------------------
755	else if (trim(modif).eq.'emis') then
756
757	print*,'new value?'
758	read(,) emisread
759
760	do i=1,ngridmx
761	emis(i)=emisread
762	enddo
763
764	c qname : change tracer name
765	c --------------------------
766	else if (trim(modif).eq.'qname') then
767	yes='y'
768	do while (yes.eq.'y')
769	write(,) 'Which tracer name do you want to change ?'
770	do iq=1,nqmx
771	write(*,'(i3,a3,a20)')iq,' : ',trim(tnom(iq))
772	enddo
773	write(*,'(a35,i3)')
774	& '(enter tracer number; between 1 and ',nqmx
775	write(,)' or any other value to quit this option)'
776	read(,) iq
777	if ((iq.ge.1).and.(iq.le.nqmx)) then
778	write(,)'Change tracer name ',trim(tnom(iq)),' to ?'
779	read(,) txt
780	tnom(iq)=txt
781	write(,)'Do you want to change another tracer name (y/n)?'
782	read(*,'(a)') yes
783	else
784	! inapropiate value of iq; quit this option
785	yes='n'
786	endif ! of if ((iq.ge.1).and.(iq.le.nqmx))
787	enddo ! of do while (yes.ne.'y')
788
789	c q=0 : set tracers to zero
790	c -------------------------
791	else if (modif(1:len_trim(modif)).eq.'q=0') then
792	c mise a 0 des q (traceurs)
793	write(,) 'Tracers set to 0 (1.E-30 in fact)'
794	DO iq =1, nqmx
795	DO l=1,llm
796	DO j=1,jjp1
797	DO i=1,iip1
798	q(i,j,l,iq)=1.e-30
799	ENDDO
800	ENDDO
801	ENDDO
802	ENDDO
803
804	c set surface tracers to zero
805	DO iq =1, nqmx
806	DO ig=1,ngridmx
807	qsurf(ig,iq)=0.
808	ENDDO
809	ENDDO
810
811	c q=x : initialise tracer manually
812	c --------------------------------
813	else if (modif(1:len_trim(modif)).eq.'q=x') then
814	write(,) 'Which tracer do you want to modify ?'
815	do iq=1,nqmx
816	write(,)iq,' : ',trim(tnom(iq))
817	enddo
818	write(,) '(choose between 1 and ',nqmx,')'
819	read(,) iq
820	write(,)'mixing ratio of tracer ',trim(tnom(iq)),
821	& ' ? (kg/kg)'
822	read(,) val
823	DO l=1,llm
824	DO j=1,jjp1
825	DO i=1,iip1
826	q(i,j,l,iq)=val
827	ENDDO
828	ENDDO
829	ENDDO
830	write(,) 'SURFACE value of tracer ',trim(tnom(iq)),
831	& ' ? (kg/m2)'
832	read(,) val
833	DO ig=1,ngridmx
834	qsurf(ig,iq)=val
835	ENDDO
836
837	c ini_q : Initialize tracers for chemistry
838	c -----------------------------------------------
839	else if (modif(1:len_trim(modif)) .eq. 'ini_q') then
840	c For more than 32 layers, possible to initiate thermosphere only
841	thermo=0
842	yes=' '
843	if (llm.gt.32) then
844	do while ((yes.ne.'y').and.(yes.ne.'n'))
845	write(,)'',
846	& 'initialisation for thermosphere only? (y/n)'
847	read(*,fmt='(a)') yes
848	if (yes.eq.'y') then
849	thermo=1
850	else
851	thermo=0
852	endif
853	enddo
854	endif
855
856	c call inichim_newstart(q,ps,sig,nqmx,latfi,lonfi,airefi,
857	c $ thermo,qsurf)
858	write(,) 'Chemical species initialized'
859
860	if (thermo.eq.0) then
861	c mise a 0 des qsurf (traceurs a la surface)
862	DO iq =1, nqmx
863	DO ig=1,ngridmx
864	qsurf(ig,iq)=0.
865	ENDDO
866	ENDDO
867	endif
868
869	c ini_q-H2O : as above except for the water vapour tracer
870	c ------------------------------------------------------
871	else if (modif(1:len_trim(modif)) .eq. 'ini_q-H2O') then
872	! for more than 32 layers, possible to initiate thermosphere only
873	thermo=0
874	yes=' '
875	if(llm.gt.32) then
876	do while ((yes.ne.'y').and.(yes.ne.'n'))
877	write(,)'',
878	& 'initialisation for thermosphere only? (y/n)'
879	read(*,fmt='(a)') yes
880	if (yes.eq.'y') then
881	thermo=1
882	else
883	thermo=0
884	endif
885	enddo
886	endif
887	c call inichim_newstart(q,ps,sig,nqmx-1,latfi,lonfi,airefi,
888	c $ thermo,qsurf)
889	c write(,) 'Initialized chem. species exept last (H2O)'
890
891	if (thermo.eq.0) then
892	c set surface tracers to zero, except water ice
893	DO iq =1, nqmx
894	if (iq.ne.igcm_h2o_ice) then
895	DO ig=1,ngridmx
896	qsurf(ig,iq)=0.
897	ENDDO
898	endif
899	ENDDO
900	endif
901
902	c ini_q-iceH2O : as above exept for ice et H2O
903	c -----------------------------------------------
904	else if (modif(1:len_trim(modif)) .eq. 'ini_q-iceH2O') then
905	c For more than 32 layers, possible to initiate thermosphere only
906	thermo=0
907	yes=' '
908	if(llm.gt.32) then
909	do while ((yes.ne.'y').and.(yes.ne.'n'))
910	write(,)'',
911	& 'initialisation for thermosphere only? (y/n)'
912	read(*,fmt='(a)') yes
913	if (yes.eq.'y') then
914	thermo=1
915	else
916	thermo=0
917	endif
918	enddo
919	endif
920
921	c call inichim_newstart(q,ps,sig,nqmx-2,latfi,lonfi,airefi,
922	c $ thermo,qsurf)
923	c write(,) 'Initialized chem. species exept ice and H2O'
924
925	if (thermo.eq.0) then
926	c set surface tracers to zero, except water ice
927	DO iq =1, nqmx
928	if (iq.ne.igcm_h2o_ice) then
929	DO ig=1,ngridmx
930	qsurf(ig,iq)=0.
931	ENDDO
932	endif
933	ENDDO
934	endif
935
936	c wetstart : wet atmosphere with a north to south gradient
937	c --------------------------------------------------------
938	else if (modif(1:len_trim(modif)) .eq. 'wetstart') then
939	! check that there is indeed a water vapor tracer
940	if (igcm_h2o_vap.eq.0) then
941	write(,) "No water vapour tracer! Can't use this option"
942	stop
943	endif
944	DO l=1,llm
945	DO j=1,jjp1
946	DO i=1,iip1
947	q(i,j,l,igcm_h2o_vap)=150.e-6 * (rlatu(j)+pi/2.) / pi
948	ENDDO
949	ENDDO
950	ENDDO
951
952	write(,) 'Water mass mixing ratio at north pole='
953	* ,q(1,1,1,igcm_h2o_vap)
954	write(,) '---------------------------south pole='
955	* ,q(1,jjp1,1,igcm_h2o_vap)
956
957	c noglacier : remove tropical water ice (to initialize high res sim)
958	c --------------------------------------------------
959	else if (modif(1:len_trim(modif)) .eq. 'noglacier') then
960	if (igcm_h2o_ice.eq.0) then
961	write(,) "No water ice tracer! Can't use this option"
962	stop
963	endif
964	do ig=1,ngridmx
965	j=(ig-2)/iim +2
966	if(ig.eq.1) j=1
967	write(,) 'OK: remove surface ice for \|lat\|<45'
968	if (abs(rlatu(j)*180./pi).lt.45.) then
969	qsurf(ig,igcm_h2o_ice)=0.
970	end if
971	end do
972
973
974	c watercapn : H20 ice on permanent northern cap
975	c --------------------------------------------------
976	else if (modif(1:len_trim(modif)) .eq. 'watercapn') then
977	if (igcm_h2o_ice.eq.0) then
978	write(,) "No water ice tracer! Can't use this option"
979	stop
980	endif
981
982	DO j=1,jjp1
983	DO i=1,iim
984	ig=1+(j-2)*iim +i
985	if(j.eq.1) ig=1
986	if(j.eq.jjp1) ig=ngridmx
987
988	if (rlatu(j)*180./pi.gt.80.) then
989	qsurf(ig,igcm_h2o_ice)=3.4e3
990	!do isoil=1,nsoilmx
991	! ith(i,j,isoil) = 18000. ! thermal inertia
992	!enddo
993	write(,)' ==> Ice mesh North boundary (deg)= ',
994	& rlatv(j-1)*180./pi
995	end if
996	ENDDO
997	ENDDO
998	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
999
1000	c$$$ do ig=1,ngridmx
1001	c$$$ j=(ig-2)/iim +2
1002	c$$$ if(ig.eq.1) j=1
1003	c$$$ if (rlatu(j)*180./pi.gt.80.) then
1004	c$$$
1005	c$$$ qsurf(ig,igcm_h2o_ice)=1.e5
1006	c$$$ qsurf(ig,igcm_h2o_vap)=0.0!1.e5
1007	c$$$
1008	c$$$ write(,) 'ig=',ig,' H2O ice mass (kg/m2)= ',
1009	c$$$ & qsurf(ig,igcm_h2o_ice)
1010	c$$$
1011	c$$$ write(,)' ==> Ice mesh South boundary (deg)= ',
1012	c$$$ & rlatv(j)*180./pi
1013	c$$$ end if
1014	c$$$ enddo
1015
1016	c watercaps : H20 ice on permanent southern cap
1017	c -------------------------------------------------
1018	else if (modif(1:len_trim(modif)) .eq. 'watercaps') then
1019	if (igcm_h2o_ice.eq.0) then
1020	write(,) "No water ice tracer! Can't use this option"
1021	stop
1022	endif
1023
1024	DO j=1,jjp1
1025	DO i=1,iim
1026	ig=1+(j-2)*iim +i
1027	if(j.eq.1) ig=1
1028	if(j.eq.jjp1) ig=ngridmx
1029
1030	if (rlatu(j)*180./pi.lt.-80.) then
1031	qsurf(ig,igcm_h2o_ice)=3.4e3
1032	!do isoil=1,nsoilmx
1033	! ith(i,j,isoil) = 18000. ! thermal inertia
1034	!enddo
1035	write(,)' ==> Ice mesh South boundary (deg)= ',
1036	& rlatv(j-1)*180./pi
1037	end if
1038	ENDDO
1039	ENDDO
1040	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
1041
1042	c$$$ do ig=1,ngridmx
1043	c$$$ j=(ig-2)/iim +2
1044	c$$$ if(ig.eq.1) j=1
1045	c$$$ if (rlatu(j)*180./pi.lt.-80.) then
1046	c$$$ qsurf(ig,igcm_h2o_ice)=1.e5
1047	c$$$ qsurf(ig,igcm_h2o_vap)=0.0 !1.e5
1048	c$$$
1049	c$$$ write(,) 'ig=',ig,' H2O ice mass (kg/m2)= ',
1050	c$$$ & qsurf(ig,igcm_h2o_ice)
1051	c$$$ write(,)' ==> Ice mesh North boundary (deg)= ',
1052	c$$$ & rlatv(j-1)*180./pi
1053	c$$$ end if
1054	c$$$ enddo
1055
1056
1057	c noacglac : H2O ice across highest terrain
1058	c --------------------------------------------
1059	else if (modif(1:len_trim(modif)) .eq. 'noacglac') then
1060	if (igcm_h2o_ice.eq.0) then
1061	write(,) "No water ice tracer! Can't use this option"
1062	stop
1063	endif
1064	DO j=1,jjp1
1065	DO i=1,iim
1066	ig=1+(j-2)*iim +i
1067	if(j.eq.1) ig=1
1068	if(j.eq.jjp1) ig=ngridmx
1069
1070	if(phis(i,j).gt.1000.*g)then
1071	alb(i,j) = 0.5 ! snow value
1072	do isoil=1,nsoilmx
1073	ith(i,j,isoil) = 18000. ! thermal inertia
1074	! this leads to rnat set to 'ocean' in physiq.F90
1075	! actually no, because it is soil not surface
1076	enddo
1077	endif
1078	ENDDO
1079	ENDDO
1080	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
1081	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
1082	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
1083
1084
1085
1086	c oborealis : H2O oceans across Vastitas Borealis
1087	c -----------------------------------------------
1088	else if (modif(1:len_trim(modif)) .eq. 'oborealis') then
1089	if (igcm_h2o_ice.eq.0) then
1090	write(,) "No water ice tracer! Can't use this option"
1091	stop
1092	endif
1093	DO j=1,jjp1
1094	DO i=1,iim
1095	ig=1+(j-2)*iim +i
1096	if(j.eq.1) ig=1
1097	if(j.eq.jjp1) ig=ngridmx
1098
1099	if(phis(i,j).lt.-4000.*g)then
1100	! if( (phis(i,j).lt.-4000.*g)
1101	! & .and. (rlatu(j)*180./pi.lt.0.) )then ! south hemisphere only
1102	! phis(i,j)=-8000.0*g ! proper ocean
1103	phis(i,j)=-4000.0*g ! scanty ocean
1104
1105	alb(i,j) = 0.07 ! oceanic value
1106	do isoil=1,nsoilmx
1107	ith(i,j,isoil) = 18000. ! thermal inertia
1108	! this leads to rnat set to 'ocean' in physiq.F90
1109	! actually no, because it is soil not surface
1110	enddo
1111	endif
1112	ENDDO
1113	ENDDO
1114	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
1115	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
1116	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
1117
1118	c iborealis : H2O ice in Northern plains
1119	c --------------------------------------
1120	else if (modif(1:len_trim(modif)) .eq. 'iborealis') then
1121	if (igcm_h2o_ice.eq.0) then
1122	write(,) "No water ice tracer! Can't use this option"
1123	stop
1124	endif
1125	DO j=1,jjp1
1126	DO i=1,iim
1127	ig=1+(j-2)*iim +i
1128	if(j.eq.1) ig=1
1129	if(j.eq.jjp1) ig=ngridmx
1130
1131	if(phis(i,j).lt.-4000.*g)then
1132	!qsurf(ig,igcm_h2o_ice)=1.e3
1133	qsurf(ig,igcm_h2o_ice)=241.4 ! to make total 33 kg m^-2
1134	endif
1135	ENDDO
1136	ENDDO
1137	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
1138	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
1139	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
1140
1141
1142	c oceanball : H2O liquid everywhere
1143	c ----------------------------
1144	else if (modif(1:len_trim(modif)) .eq. 'oceanball') then
1145	if (igcm_h2o_ice.eq.0) then
1146	write(,) "No water ice tracer! Can't use this option"
1147	stop
1148	endif
1149	DO j=1,jjp1
1150	DO i=1,iim
1151	ig=1+(j-2)*iim +i
1152	if(j.eq.1) ig=1
1153	if(j.eq.jjp1) ig=ngridmx
1154
1155	qsurf(ig,igcm_h2o_vap)=0.0 ! for ocean, this is infinite source
1156	qsurf(ig,igcm_h2o_ice)=0.0
1157	alb(i,j) = 0.07 ! ocean value
1158
1159	do isoil=1,nsoilmx
1160	ith(i,j,isoil) = 18000. ! thermal inertia
1161	!ith(i,j,isoil) = 50. ! extremely low for test
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	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
1170
1171	c iceball : H2O ice everywhere
1172	c ----------------------------
1173	else if (modif(1:len_trim(modif)) .eq. 'iceball') then
1174	if (igcm_h2o_ice.eq.0) then
1175	write(,) "No water ice tracer! Can't use this option"
1176	stop
1177	endif
1178	DO j=1,jjp1
1179	DO i=1,iim
1180	ig=1+(j-2)*iim +i
1181	if(j.eq.1) ig=1
1182	if(j.eq.jjp1) ig=ngridmx
1183
1184	qsurf(ig,igcm_h2o_vap)=-50. ! for ocean, this is infinite source
1185	qsurf(ig,igcm_h2o_ice)=50. ! == to 0.05 m of oceanic ice
1186	alb(i,j) = 0.6 ! ice albedo value
1187
1188	do isoil=1,nsoilmx
1189	!ith(i,j,isoil) = 18000. ! thermal inertia
1190	! this leads to rnat set to 'ocean' in physiq.F90
1191	enddo
1192
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
1198	c isotherm : Isothermal temperatures and no winds
1199	c -----------------------------------------------
1200	else if (modif(1:len_trim(modif)) .eq. 'isotherm') then
1201
1202	write(,)'Isothermal temperature of the atmosphere,
1203	& surface and subsurface'
1204	write(,) 'Value of this temperature ? :'
1205	203 read(,,iostat=ierr) Tiso
1206	if(ierr.ne.0) goto 203
1207
1208	do ig=1, ngridmx
1209	tsurf(ig) = Tiso
1210	end do
1211	do l=2,nsoilmx
1212	do ig=1, ngridmx
1213	tsoil(ig,l) = Tiso
1214	end do
1215	end do
1216	DO j=1,jjp1
1217	DO i=1,iim
1218	Do l=1,llm
1219	Tset(i,j,l)=Tiso
1220	end do
1221	end do
1222	end do
1223	flagtset=.true.
1224	call initial0(llm*ip1jmp1,ucov)
1225	call initial0(llm*ip1jm,vcov)
1226	call initial0(ngridmx*(llm+1),q2)
1227
1228	c radequi : Radiative equilibrium profile of temperatures and no winds
1229	c --------------------------------------------------------------------
1230	else if (modif(1:len_trim(modif)) .eq. 'radequi') then
1231
1232	write(,)'radiative equilibrium temperature profile'
1233
1234	do ig=1, ngridmx
1235	teque= 335.0-60.0sin(latfi(ig))sin(latfi(ig))-
1236	& 10.0cos(latfi(ig))cos(latfi(ig))
1237	tsurf(ig) = MAX(220.0,teque)
1238	end do
1239	do l=2,nsoilmx
1240	do ig=1, ngridmx
1241	tsoil(ig,l) = tsurf(ig)
1242	end do
1243	end do
1244	DO j=1,jjp1
1245	DO i=1,iim
1246	Do l=1,llm
1247	teque=335.-60.0sin(rlatu(j))sin(rlatu(j))-
1248	& 10.0cos(rlatu(j))cos(rlatu(j))
1249	Tset(i,j,l)=MAX(220.0,teque)
1250	end do
1251	end do
1252	end do
1253	flagtset=.true.
1254	call initial0(llm*ip1jmp1,ucov)
1255	call initial0(llm*ip1jm,vcov)
1256	call initial0(ngridmx*(llm+1),q2)
1257
1258	c coldstart : T set 1K above CO2 frost point and no winds
1259	c ------------------------------------------------
1260	else if (modif(1:len_trim(modif)) .eq. 'coldstart') then
1261
1262	write(,)'set temperature of the atmosphere,'
1263	&,'surface and subsurface how many degrees above CO2 frost point?'
1264	204 read(,,iostat=ierr) Tabove
1265	if(ierr.ne.0) goto 204
1266
1267	DO j=1,jjp1
1268	DO i=1,iim
1269	ig=1+(j-2)*iim +i
1270	if(j.eq.1) ig=1
1271	if(j.eq.jjp1) ig=ngridmx
1272	tsurf(ig) = (-3167.8)/(log(.01*ps(i,j))-23.23)+Tabove
1273	END DO
1274	END DO
1275	do l=1,nsoilmx
1276	do ig=1, ngridmx
1277	tsoil(ig,l) = tsurf(ig)
1278	end do
1279	end do
1280	DO j=1,jjp1
1281	DO i=1,iim
1282	Do l=1,llm
1283	pp = aps(l) +bps(l)*ps(i,j)
1284	Tset(i,j,l)=(-3167.8)/(log(.01*pp)-23.23)+Tabove
1285	end do
1286	end do
1287	end do
1288
1289	flagtset=.true.
1290	call initial0(llm*ip1jmp1,ucov)
1291	call initial0(llm*ip1jm,vcov)
1292	call initial0(ngridmx*(llm+1),q2)
1293
1294
1295	c co2ice=0 : remove CO2 polar ice caps'
1296	c ------------------------------------------------
1297	else if (modif(1:len_trim(modif)) .eq. 'co2ice=0') then
1298	! check that there is indeed a co2_ice tracer ...
1299	if (igcm_co2_ice.ne.0) then
1300	do ig=1,ngridmx
1301	!co2ice(ig)=0
1302	qsurf(ig,igcm_co2_ice)=0
1303	emis(ig)=emis(ngridmx/2)
1304	end do
1305	else
1306	write(,) "Can't remove CO2 ice!! (no co2_ice tracer)"
1307	endif
1308
1309	! therm_ini_s: (re)-set soil thermal inertia to reference surface values
1310	! ----------------------------------------------------------------------
1311
1312	else if (modif(1:len_trim(modif)).eq.'therm_ini_s') then
1313	! write(,)"surfithfi(1):",surfithfi(1)
1314	do isoil=1,nsoilmx
1315	inertiedat(1:ngridmx,isoil)=surfithfi(1:ngridmx)
1316	enddo
1317	write(,)'OK: Soil thermal inertia has been reset to referenc
1318	&e surface values'
1319	! write(,)"inertiedat(1,1):",inertiedat(1,1)
1320	ithfi(:,:)=inertiedat(:,:)
1321	! recast ithfi() onto ith()
1322	call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
1323	! Check:
1324	! do i=1,iip1
1325	! do j=1,jjp1
1326	! do isoil=1,nsoilmx
1327	! write(77,*) i,j,isoil," ",ith(i,j,isoil)
1328	! enddo
1329	! enddo
1330	! enddo
1331
1332
1333
1334
1335	c$$$! subsoilice_n: Put deep ice layer in northern hemisphere soil
1336	c$$$! ------------------------------------------------------------
1337	c$$$
1338	c$$$ else if (modif(1:len_trim(modif)).eq.'subsoilice_n') then
1339	c$$$
1340	c$$$ write(,)'From which latitude (in deg.), up to the north pole,
1341	c$$$ &should we put subterranean ice?'
1342	c$$$ ierr=1
1343	c$$$ do while (ierr.ne.0)
1344	c$$$ read(,,iostat=ierr) val
1345	c$$$ if (ierr.eq.0) then ! got a value
1346	c$$$ ! do a sanity check
1347	c$$$ if((val.lt.0.).or.(val.gt.90)) then
1348	c$$$ write(,)'Latitude should be between 0 and 90 deg. !!!'
1349	c$$$ ierr=1
1350	c$$$ else ! find corresponding jref (nearest latitude)
1351	c$$$ ! note: rlatu(:) contains decreasing values of latitude
1352	c$$$ ! starting from PI/2 to -PI/2
1353	c$$$ do j=1,jjp1
1354	c$$$ if ((rlatu(j)*180./pi.ge.val).and.
1355	c$$$ & (rlatu(j+1)*180./pi.le.val)) then
1356	c$$$ ! find which grid point is nearest to val:
1357	c$$$ if (abs(rlatu(j)*180./pi-val).le.
1358	c$$$ & abs((rlatu(j+1)*180./pi-val))) then
1359	c$$$ jref=j
1360	c$$$ else
1361	c$$$ jref=j+1
1362	c$$$ endif
1363	c$$$
1364	c$$$ write(,)'Will use nearest grid latitude which is:',
1365	c$$$ & rlatu(jref)*180./pi
1366	c$$$ endif
1367	c$$$ enddo ! of do j=1,jjp1
1368	c$$$ endif ! of if((val.lt.0.).or.(val.gt.90))
1369	c$$$ endif !of if (ierr.eq.0)
1370	c$$$ enddo ! of do while
1371	c$$$
1372	c$$$ ! Build layers() (as in soil_settings.F)
1373	c$$$ val2=sqrt(mlayer(0)*mlayer(1))
1374	c$$$ val3=mlayer(1)/mlayer(0)
1375	c$$$ do isoil=1,nsoilmx
1376	c$$$ layer(isoil)=val2(val3*(isoil-1))
1377	c$$$ enddo
1378	c$$$
1379	c$$$ write(,)'At which depth (in m.) does the ice layer begin?'
1380	c$$$ write(,)'(currently, the deepest soil layer extends down to:'
1381	c$$$ & ,layer(nsoilmx),')'
1382	c$$$ ierr=1
1383	c$$$ do while (ierr.ne.0)
1384	c$$$ read(,,iostat=ierr) val2
1385	c$$$! write(,)'val2:',val2,'ierr=',ierr
1386	c$$$ if (ierr.eq.0) then ! got a value, but do a sanity check
1387	c$$$ if(val2.gt.layer(nsoilmx)) then
1388	c$$$ write(,)'Depth should be less than ',layer(nsoilmx)
1389	c$$$ ierr=1
1390	c$$$ endif
1391	c$$$ if(val2.lt.layer(1)) then
1392	c$$$ write(,)'Depth should be more than ',layer(1)
1393	c$$$ ierr=1
1394	c$$$ endif
1395	c$$$ endif
1396	c$$$ enddo ! of do while
1397	c$$$
1398	c$$$ ! find the reference index iref the depth corresponds to
1399	c$$$! if (val2.lt.layer(1)) then
1400	c$$$! iref=1
1401	c$$$! else
1402	c$$$ do isoil=1,nsoilmx-1
1403	c$$$ if((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1)))
1404	c$$$ & then
1405	c$$$ iref=isoil
1406	c$$$ exit
1407	c$$$ endif
1408	c$$$ enddo
1409	c$$$! endif
1410	c$$$
1411	c$$$! write(,)'iref:',iref,' jref:',jref
1412	c$$$! write(,)'layer',layer
1413	c$$$! write(,)'mlayer',mlayer
1414	c$$$
1415	c$$$ ! thermal inertia of the ice:
1416	c$$$ ierr=1
1417	c$$$ do while (ierr.ne.0)
1418	c$$$ write(,)'What is the value of subterranean ice thermal inert
1419	c$$$ &ia? (e.g.: 2000)'
1420	c$$$ read(,,iostat=ierr)iceith
1421	c$$$ enddo ! of do while
1422	c$$$
1423	c$$$ ! recast ithfi() onto ith()
1424	c$$$ call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
1425	c$$$
1426	c$$$ do j=1,jref
1427	c$$$! write(,)'j:',j,'rlatu(j)180./pi:',rlatu(j)180./pi
1428	c$$$ do i=1,iip1 ! loop on longitudes
1429	c$$$ ! Build "equivalent" thermal inertia for the mixed layer
1430	c$$$ ith(i,j,iref+1)=sqrt((layer(iref+1)-layer(iref))/
1431	c$$$ & (((val2-layer(iref))/(ith(i,j,iref)**2))+
1432	c$$$ & ((layer(iref+1)-val2)/(iceith)**2)))
1433	c$$$ ! Set thermal inertia of lower layers
1434	c$$$ do isoil=iref+2,nsoilmx
1435	c$$$ ith(i,j,isoil)=iceith ! ice
1436	c$$$ enddo
1437	c$$$ enddo ! of do i=1,iip1
1438	c$$$ enddo ! of do j=1,jjp1
1439	c$$$
1440	c$$$
1441	c$$$ CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
1442	c$$$
1443	c$$$! do i=1,nsoilmx
1444	c$$$! write(,)'i:',i,'ithfi(1,i):',ithfi(1,i)
1445	c$$$! enddo
1446
1447
1448	c$$$! subsoilice_s: Put deep ice layer in southern hemisphere soil
1449	c$$$! ------------------------------------------------------------
1450	c$$$
1451	c$$$ else if (modif(1:len_trim(modif)).eq.'subsoilice_s') then
1452	c$$$
1453	c$$$ write(,)'From which latitude (in deg.), down to the south pol
1454	c$$$ &e, should we put subterranean ice?'
1455	c$$$ ierr=1
1456	c$$$ do while (ierr.ne.0)
1457	c$$$ read(,,iostat=ierr) val
1458	c$$$ if (ierr.eq.0) then ! got a value
1459	c$$$ ! do a sanity check
1460	c$$$ if((val.gt.0.).or.(val.lt.-90)) then
1461	c$$$ write(,)'Latitude should be between 0 and -90 deg. !!!'
1462	c$$$ ierr=1
1463	c$$$ else ! find corresponding jref (nearest latitude)
1464	c$$$ ! note: rlatu(:) contains decreasing values of latitude
1465	c$$$ ! starting from PI/2 to -PI/2
1466	c$$$ do j=1,jjp1
1467	c$$$ if ((rlatu(j)*180./pi.ge.val).and.
1468	c$$$ & (rlatu(j+1)*180./pi.le.val)) then
1469	c$$$ ! find which grid point is nearest to val:
1470	c$$$ if (abs(rlatu(j)*180./pi-val).le.
1471	c$$$ & abs((rlatu(j+1)*180./pi-val))) then
1472	c$$$ jref=j
1473	c$$$ else
1474	c$$$ jref=j+1
1475	c$$$ endif
1476	c$$$
1477	c$$$ write(,)'Will use nearest grid latitude which is:',
1478	c$$$ & rlatu(jref)*180./pi
1479	c$$$ endif
1480	c$$$ enddo ! of do j=1,jjp1
1481	c$$$ endif ! of if((val.lt.0.).or.(val.gt.90))
1482	c$$$ endif !of if (ierr.eq.0)
1483	c$$$ enddo ! of do while
1484	c$$$
1485	c$$$ ! Build layers() (as in soil_settings.F)
1486	c$$$ val2=sqrt(mlayer(0)*mlayer(1))
1487	c$$$ val3=mlayer(1)/mlayer(0)
1488	c$$$ do isoil=1,nsoilmx
1489	c$$$ layer(isoil)=val2(val3*(isoil-1))
1490	c$$$ enddo
1491	c$$$
1492	c$$$ write(,)'At which depth (in m.) does the ice layer begin?'
1493	c$$$ write(,)'(currently, the deepest soil layer extends down to:'
1494	c$$$ & ,layer(nsoilmx),')'
1495	c$$$ ierr=1
1496	c$$$ do while (ierr.ne.0)
1497	c$$$ read(,,iostat=ierr) val2
1498	c$$$! write(,)'val2:',val2,'ierr=',ierr
1499	c$$$ if (ierr.eq.0) then ! got a value, but do a sanity check
1500	c$$$ if(val2.gt.layer(nsoilmx)) then
1501	c$$$ write(,)'Depth should be less than ',layer(nsoilmx)
1502	c$$$ ierr=1
1503	c$$$ endif
1504	c$$$ if(val2.lt.layer(1)) then
1505	c$$$ write(,)'Depth should be more than ',layer(1)
1506	c$$$ ierr=1
1507	c$$$ endif
1508	c$$$ endif
1509	c$$$ enddo ! of do while
1510	c$$$
1511	c$$$ ! find the reference index iref the depth corresponds to
1512	c$$$ do isoil=1,nsoilmx-1
1513	c$$$ if((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1)))
1514	c$$$ & then
1515	c$$$ iref=isoil
1516	c$$$ exit
1517	c$$$ endif
1518	c$$$ enddo
1519	c$$$
1520	c$$$! write(,)'iref:',iref,' jref:',jref
1521	c$$$
1522	c$$$ ! thermal inertia of the ice:
1523	c$$$ ierr=1
1524	c$$$ do while (ierr.ne.0)
1525	c$$$ write(,)'What is the value of subterranean ice thermal inert
1526	c$$$ &ia? (e.g.: 2000)'
1527	c$$$ read(,,iostat=ierr)iceith
1528	c$$$ enddo ! of do while
1529	c$$$
1530	c$$$ ! recast ithfi() onto ith()
1531	c$$$ call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
1532	c$$$
1533	c$$$ do j=jref,jjp1
1534	c$$$! write(,)'j:',j,'rlatu(j)180./pi:',rlatu(j)180./pi
1535	c$$$ do i=1,iip1 ! loop on longitudes
1536	c$$$ ! Build "equivalent" thermal inertia for the mixed layer
1537	c$$$ ith(i,j,iref+1)=sqrt((layer(iref+1)-layer(iref))/
1538	c$$$ & (((val2-layer(iref))/(ith(i,j,iref)**2))+
1539	c$$$ & ((layer(iref+1)-val2)/(iceith)**2)))
1540	c$$$ ! Set thermal inertia of lower layers
1541	c$$$ do isoil=iref+2,nsoilmx
1542	c$$$ ith(i,j,isoil)=iceith ! ice
1543	c$$$ enddo
1544	c$$$ enddo ! of do i=1,iip1
1545	c$$$ enddo ! of do j=jref,jjp1
1546	c$$$
1547	c$$$
1548	c$$$ CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
1549
1550
1551	c$$$c 'mons_ice' : use MONS data to build subsurface ice table
1552	c$$$c --------------------------------------------------------
1553	c$$$ else if (modif(1:len_trim(modif)).eq.'mons_ice') then
1554	c$$$
1555	c$$$ ! 1. Load MONS data
1556	c$$$ call load_MONS_data(MONS_Hdn,MONS_d21)
1557	c$$$
1558	c$$$ ! 2. Get parameters from user
1559	c$$$ ierr=1
1560	c$$$ do while (ierr.ne.0)
1561	c$$$ write(,) "Coefficient to apply to MONS 'depth' in Northern",
1562	c$$$ & " Hemisphere?"
1563	c$$$ write(,) " (should be somewhere between 3.2e-4 and 1.3e-3)"
1564	c$$$ read(,,iostat=ierr) MONS_coeffN
1565	c$$$ enddo
1566	c$$$ ierr=1
1567	c$$$ do while (ierr.ne.0)
1568	c$$$ write(,) "Coefficient to apply to MONS 'depth' in Southern",
1569	c$$$ & " Hemisphere?"
1570	c$$$ write(,) " (should be somewhere between 3.2e-4 and 1.3e-3)"
1571	c$$$ read(,,iostat=ierr) MONS_coeffS
1572	c$$$ enddo
1573	c$$$ ierr=1
1574	c$$$ do while (ierr.ne.0)
1575	c$$$ write(,) "Value of subterranean ice thermal inertia?"
1576	c$$$ write(,) " (e.g.: 2000, or perhaps 2290)"
1577	c$$$ read(,,iostat=ierr) iceith
1578	c$$$ enddo
1579	c$$$
1580	c$$$ ! 3. Build subterranean thermal inertia
1581	c$$$
1582	c$$$ ! initialise subsurface inertia with reference surface values
1583	c$$$ do isoil=1,nsoilmx
1584	c$$$ ithfi(1:ngridmx,isoil)=surfithfi(1:ngridmx)
1585	c$$$ enddo
1586	c$$$ ! recast ithfi() onto ith()
1587	c$$$ call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
1588	c$$$
1589	c$$$ do i=1,iip1 ! loop on longitudes
1590	c$$$ do j=1,jjp1 ! loop on latitudes
1591	c$$$ ! set MONS_coeff
1592	c$$$ if (rlatu(j).ge.0) then ! northern hemisphere
1593	c$$$ ! N.B: rlatu(:) contains decreasing values of latitude
1594	c$$$ ! starting from PI/2 to -PI/2
1595	c$$$ MONS_coeff=MONS_coeffN
1596	c$$$ else ! southern hemisphere
1597	c$$$ MONS_coeff=MONS_coeffS
1598	c$$$ endif
1599	c$$$ ! check if we should put subterranean ice
1600	c$$$ if (MONS_Hdn(i,j).ge.14.0) then ! no ice if Hdn<14%
1601	c$$$ ! compute depth at which ice lies:
1602	c$$$ val=MONS_d21(i,j)*MONS_coeff
1603	c$$$ ! compute val2= the diurnal skin depth of surface inertia
1604	c$$$ ! assuming a volumetric heat cap. of C=1.e6 J.m-3.K-1
1605	c$$$ val2=ith(i,j,1)1.e-6sqrt(88775./3.14159)
1606	c$$$ if (val.lt.val2) then
1607	c$$$ ! ice must be below the (surface inertia) diurnal skin depth
1608	c$$$ val=val2
1609	c$$$ endif
1610	c$$$ if (val.lt.layer(nsoilmx)) then ! subterranean ice
1611	c$$$ ! find the reference index iref that depth corresponds to
1612	c$$$ iref=0
1613	c$$$ do isoil=1,nsoilmx-1
1614	c$$$ if ((val.ge.layer(isoil)).and.(val.lt.layer(isoil+1)))
1615	c$$$ & then
1616	c$$$ iref=isoil
1617	c$$$ exit
1618	c$$$ endif
1619	c$$$ enddo
1620	c$$$ ! Build "equivalent" thermal inertia for the mixed layer
1621	c$$$ ith(i,j,iref+1)=sqrt((layer(iref+1)-layer(iref))/
1622	c$$$ & (((val-layer(iref))/(ith(i,j,iref+1)**2))+
1623	c$$$ & ((layer(iref+1)-val)/(iceith)**2)))
1624	c$$$ ! Set thermal inertia of lower layers
1625	c$$$ do isoil=iref+2,nsoilmx
1626	c$$$ ith(i,j,isoil)=iceith
1627	c$$$ enddo
1628	c$$$ endif ! of if (val.lt.layer(nsoilmx))
1629	c$$$ endif ! of if (MONS_Hdn(i,j).lt.14.0)
1630	c$$$ enddo ! do j=1,jjp1
1631	c$$$ enddo ! do i=1,iip1
1632	c$$$
1633	c$$$! Check:
1634	c$$$! do i=1,iip1
1635	c$$$! do j=1,jjp1
1636	c$$$! do isoil=1,nsoilmx
1637	c$$$! write(77,*) i,j,isoil," ",ith(i,j,isoil)
1638	c$$$! enddo
1639	c$$$! enddo
1640	c$$$! enddo
1641	c$$$
1642	c$$$ ! recast ith() into ithfi()
1643	c$$$ CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
1644	c$$$
1645	c$$$ else
1646	c$$$ write(,) ' Unknown (misspelled?) option!!!'
1647	end if ! of if (modif(1:len_trim(modif)) .eq. '...') elseif ...
1648
1649
1650
1651
1652
1653
1654
1655
1656
1657
1658
1659	enddo ! of do ! infinite loop on liste of changes
1660
1661	999 continue
1662
1663
1664	c=======================================================================
1665	c Initialisation for cloud fraction and oceanic ice (added by BC 2010)
1666	c=======================================================================
1667	DO ig=1, ngridmx
1668	totalfrac(ig)=0.5
1669	DO l=1,llm
1670	cloudfrac(ig,l)=0.5
1671	ENDDO
1672	ENDDO
1673
1674	c========================================================
1675
1676	! DO ig=1,ngridmx
1677	! IF(tsurf(ig) .LT. 273.16-1.8) THEN
1678	! hice(ig)=(273.16-1.8-tsurf(ig))/(273.16-1.8-240)*1.
1679	! hice(ig)=min(hice(ig),1.0)
1680	! ENDIF
1681	! ENDDO
1682
1683
1684	c=======================================================================
1685	c Correct pressure on the new grid (menu 0)
1686	c=======================================================================
1687
1688	if ((choix_1.eq.0).and.(.not.flagps0)) then
1689	r = 1000.*8.31/mugaz
1690
1691	phi0=0.0
1692	do j=1,jjp1
1693	do i=1,iip1
1694	phi0 = phi0+phis(i,j)*aire(i,j)
1695	end do
1696	end do
1697	phi0=phi0/airetot
1698
1699	do j=1,jjp1
1700	do i=1,iip1
1701	ps(i,j) = ps(i,j) *
1702	! . exp((phisold_newgrid(i,j)-phis(i,j)) /
1703	. exp((phi0-phis(i,j)) /
1704	. (t(i,j,1) * r))
1705	end do
1706	end do
1707
1708	! we must renormalise pressures AGAIN, because exp(-phi) is nonlinear
1709	ptot=0.0
1710	do j=1,jjp1
1711	do i=1,iip1
1712	ptot=ptot+ps(i,j)*aire(i,j)
1713	enddo
1714	enddo
1715	do j=1,jjp1
1716	do i=1,iip1
1717	ps(i,j)=ps(i,j)*ptotn/ptot
1718	enddo
1719	enddo
1720
1721	! periodicity of surface ps in longitude
1722	do j=1,jjp1
1723	ps(1,j) = ps(iip1,j)
1724	end do
1725
1726	end if
1727
1728
1729	c=======================================================================
1730	c=======================================================================
1731
1732	c=======================================================================
1733	c Initialisation de la physique / ecriture de newstartfi :
1734	c=======================================================================
1735
1736	CALL inifilr
1737	CALL pression(ip1jmp1, ap, bp, ps, p3d)
1738
1739	c-----------------------------------------------------------------------
1740	c Initialisation pks:
1741	c-----------------------------------------------------------------------
1742
1743	CALL exner_hyb(ip1jmp1, ps, p3d, beta, pks, pk, pkf)
1744	! Calcul de la temperature potentielle teta
1745
1746	if (flagtset) then
1747	DO l=1,llm
1748	DO j=1,jjp1
1749	DO i=1,iim
1750	teta(i,j,l) = Tset(i,j,l) * cpp/pk(i,j,l)
1751	ENDDO
1752	teta (iip1,j,l)= teta (1,j,l)
1753	ENDDO
1754	ENDDO
1755	else if (choix_1.eq.0) then
1756	DO l=1,llm
1757	DO j=1,jjp1
1758	DO i=1,iim
1759	teta(i,j,l) = t(i,j,l) * cpp/pk(i,j,l)
1760	ENDDO
1761	teta (iip1,j,l)= teta (1,j,l)
1762	ENDDO
1763	ENDDO
1764	endif
1765
1766
1767	C Calcul intermediaire
1768	c
1769	if (choix_1.eq.0) then
1770	CALL massdair( p3d, masse )
1771	c
1772	print *,' ALPHAX ',alphax
1773
1774	DO l = 1, llm
1775	DO i = 1, iim
1776	xppn(i) = aire( i, 1 ) * masse( i , 1 , l )
1777	xpps(i) = aire( i,jjp1 ) * masse( i , jjp1 , l )
1778	ENDDO
1779	xpn = SUM(xppn)/apoln
1780	xps = SUM(xpps)/apols
1781	DO i = 1, iip1
1782	masse( i , 1 , l ) = xpn
1783	masse( i , jjp1 , l ) = xps
1784	ENDDO
1785	ENDDO
1786	endif
1787	phis(iip1,:) = phis(1,:)
1788
1789	CALL inidissip ( lstardis, nitergdiv, nitergrot, niterh,
1790	* tetagdiv, tetagrot , tetatemp )
1791	itau=0
1792	if (choix_1.eq.0) then
1793	day_ini=int(date)
1794	endif
1795	c
1796	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
1797
1798	CALL caldyn0( itau,ucov,vcov,teta,ps,masse,pk,phis ,
1799	* phi,w, pbaru,pbarv,day_ini+time )
1800	c CALL caldyn
1801	c $ ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
1802	c $ phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, day_ini )
1803
1804	CALL dynredem0("restart.nc",day_ini,anneeref,phis,nqmx)
1805	CALL dynredem1("restart.nc",0.0,vcov,ucov,teta,q,nqmx,masse,ps)
1806	C
1807	C Ecriture etat initial physique
1808	C
1809
1810	! do ig=1,ngridmx
1811	! print*,'surface ice in newstart=',qsurf(ig,igcm_h2o_ice)
1812	! print*,'surface liq in newstart=',qsurf(ig,igcm_h2o_vap)
1813	! enddo
1814	! stop
1815
1816
1817	call physdem1("restartfi.nc",lonfi,latfi,nsoilmx,nqmx,
1818	. dtphys,float(day_ini),
1819	. time,tsurf,tsoil,emis,q2,qsurf,
1820	. airefi,albfi,ithfi,zmea,zstd,zsig,zgam,zthe,
1821	. cloudfrac,totalfrac,hice)
1822
1823	c=======================================================================
1824	c Formats
1825	c=======================================================================
1826
1827	1 FORMAT(//10x,'la valeur de im =',i4,2x,'lue sur le fichier de dema
1828	*rrage est differente de la valeur parametree iim =',i4//)
1829	2 FORMAT(//10x,'la valeur de jm =',i4,2x,'lue sur le fichier de dema
1830	*rrage est differente de la valeur parametree jjm =',i4//)
1831	3 FORMAT(//10x,'la valeur de lllm =',i4,2x,'lue sur le fichier demar
1832	*rage est differente de la valeur parametree llm =',i4//)
1833
1834	end
1835
1836	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1837	subroutine load_MONS_data(MONS_Hdn,MONS_d21)
1838	use datafile_mod, only: datadir
1839	implicit none
1840	! routine to load Benedicte Diez MONS dataset, fill in date in southern
1841	! polar region, and interpolate the result onto the GCM grid
1842	#include"dimensions.h"
1843	#include"paramet.h"
1844	#include"comgeom.h"
1845	! arguments:
1846	real,intent(out) :: MONS_Hdn(iip1,jjp1) ! Hdn: %WEH=Mass fraction of H2O
1847	real,intent(out) :: MONS_d21(iip1,jjp1) ! ice table "depth" (in kg/m2)
1848	! N.B MONS datasets should be of dimension (iip1,jjp1)
1849	! local variables:
1850	character(len=88) :: filename="results_MONS_lat_lon_H_depth.txt"
1851	character(len=88) :: txt ! to store some text
1852	integer :: ierr,i,j
1853	integer,parameter :: nblon=180 ! number of longitudes of MONS datasets
1854	integer,parameter :: nblat=90 ! number of latitudes of MONS datasets
1855	real :: pi
1856	real :: longitudes(nblon) ! MONS dataset longitudes
1857	real :: latitudes(nblat) ! MONS dataset latitudes
1858	! MONS dataset: mass fraction of H2O where H is assumed to be in H2O
1859	real :: Hdn(nblon,nblat)
1860	real :: d21(nblon,nblat)! MONS dataset "depth" (g/cm2)
1861
1862	! Extended MONS dataset (for interp_horiz)
1863	real :: Hdnx(nblon+1,nblat)
1864	real :: d21x(nblon+1,nblat)
1865	real :: lon_bound(nblon+1) ! longitude boundaries
1866	real :: lat_bound(nblat-1) ! latitude boundaries
1867
1868	! 1. Initializations:
1869
1870	write(,) "Loading MONS data"
1871
1872	! Open MONS datafile:
1873	open(42,file=trim(datadir)//"/"//trim(filename),
1874	& status="old",iostat=ierr)
1875	if (ierr/=0) then
1876	write(,) "Error in load_MONS_data:"
1877	write(,) "Failed opening file ",
1878	& trim(datadir)//"/"//trim(filename)
1879	write(,)'Check that your path to datagcm:',trim(datadir)
1880	write(,)' is correct. You can change it in callphys.def with:'
1881	write(,)' datadir = /absolute/path/to/datagcm'
1882	write(,)'If necessary ',trim(filename),
1883	& ' (and other datafiles)'
1884	write(,)' can be obtained online at:'
1885	write(,)' http://www.lmd.jussieu.fr/~forget/datagcm/datafile'
1886	CALL ABORT
1887	else ! skip first line of file (dummy read)
1888	read(42,*) txt
1889	endif
1890
1891	pi=2.*asin(1.)
1892
1893	!2. Load MONS data (on MONS grid)
1894	do j=1,nblat
1895	do i=1,nblon
1896	! swap latitude index so latitudes go from north pole to south pole:
1897	read(42,*) latitudes(nblat-j+1),longitudes(i),
1898	& Hdn(i,nblat-j+1),d21(i,nblat-j+1)
1899	! multiply d21 by 10 to convert from g/cm2 to kg/m2
1900	d21(i,nblat-j+1)=d21(i,nblat-j+1)*10.0
1901	enddo
1902	enddo
1903	close(42)
1904
1905	! there is unfortunately no d21 data for latitudes -77 to -90
1906	! so we build some by linear interpolation between values at -75
1907	! and assuming d21=0 at the pole
1908	do j=84,90 ! latitudes(84)=-77 ; latitudes(83)=-75
1909	do i=1,nblon
1910	d21(i,j)=d21(i,83)*((latitudes(j)+90)/15.0)
1911	enddo
1912	enddo
1913
1914	! 3. Build extended MONS dataset & boundaries (for interp_horiz)
1915	! longitude boundaries (in radians):
1916	do i=1,nblon
1917	! NB: MONS data is every 2 degrees in longitude
1918	lon_bound(i)=(longitudes(i)+1.0)*pi/180.0
1919	enddo
1920	! extra 'modulo' value
1921	lon_bound(nblon+1)=lon_bound(1)+2.0*pi
1922
1923	! latitude boundaries (in radians):
1924	do j=1,nblat-1
1925	! NB: Mons data is every 2 degrees in latitude
1926	lat_bound(j)=(latitudes(j)-1.0)*pi/180.0
1927	enddo
1928
1929	! MONS datasets:
1930	do j=1,nblat
1931	Hdnx(1:nblon,j)=Hdn(1:nblon,j)
1932	Hdnx(nblon+1,j)=Hdnx(1,j)
1933	d21x(1:nblon,j)=d21(1:nblon,j)
1934	d21x(nblon+1,j)=d21x(1,j)
1935	enddo
1936
1937	! Interpolate onto GCM grid
1938	call interp_horiz(Hdnx,MONS_Hdn,nblon,nblat-1,iim,jjm,1,
1939	& lon_bound,lat_bound,rlonu,rlatv)
1940	call interp_horiz(d21x,MONS_d21,nblon,nblat-1,iim,jjm,1,
1941	& lon_bound,lat_bound,rlonu,rlatv)
1942
1943	end subroutine

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

Download in other formats:

Original Format