Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

newstart.F @ 86

Last change on this file since 86 was 38, checked in by emillour, 14 years ago
Ajout du modè Martien (mon LMDZ.MARS.BETA, du 28/01/2011) dans le rértoire mars, pour pouvoir suivre plus facilement les modifs. EM
File size: 50.8 KB

Rev	Line
[38]	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 "dimradmars.h"
	24	#include "yomaer.h"
	25	#include "planete.h"
	26	#include "paramet.h"
	27	#include "comconst.h"
	28	#include "comvert.h"
	29	#include "comgeom2.h"
	30	#include "control.h"
	31	#include "logic.h"
	32	#include "description.h"
	33	#include "ener.h"
	34	#include "temps.h"
	35	#include "lmdstd.h"
	36	#include "comdissnew.h"
	37	#include "clesph0.h"
	38	#include "serre.h"
	39	#include "netcdf.inc"
	40	#include"advtrac.h"
	41	#include"tracer.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	c et autres:
	51	c----------
	52	INTEGER lnblnk
	53	EXTERNAL lnblnk
	54
	55	c Variables pour les lectures NetCDF des fichiers "start_archive"
	56	c--------------------------------------------------
	57	INTEGER nid_dyn, nid_fi,nid,nvarid
	58	! INTEGER size
	59	INTEGER length
	60	parameter (length = 100)
	61	INTEGER tab0
	62	INTEGER NB_ETATMAX
	63	parameter (NB_ETATMAX = 100)
	64
	65	REAL date
	66	REAL p_rad,p_omeg,p_g,p_mugaz,p_daysec
	67
	68	c Variable histoire
	69	c------------------
	70	REAL vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm) ! vents covariants
	71	REAL phis(iip1,jjp1)
	72	REAL q(iip1,jjp1,llm,nqmx) ! champs advectes
	73
	74	c autre variables dynamique nouvelle grille
	75	c------------------------------------------
	76	REAL pks(iip1,jjp1)
	77	REAL w(iip1,jjp1,llm+1)
	78	REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
	79	! REAL dv(ip1jm,llm),du(ip1jmp1,llm)
	80	! REAL dh(ip1jmp1,llm),dp(ip1jmp1)
	81	REAL phi(iip1,jjp1,llm)
	82
	83	integer klatdat,klongdat
	84	PARAMETER (klatdat=180,klongdat=360)
	85
	86	c Physique sur grille scalaire
	87	c----------------------------
	88	real zmeaS(iip1,jjp1),zstdS(iip1,jjp1)
	89	real zsigS(iip1,jjp1),zgamS(iip1,jjp1),ztheS(iip1,jjp1)
	90
	91	c variable physique
	92	c------------------
	93	REAL tsurf(ngridmx) ! surface temperature
	94	REAL tsoil(ngridmx,nsoilmx) ! soil temperature
	95	REAL co2ice(ngridmx) ! CO2 ice layer
	96	REAL emis(ngridmx) ! surface emissivity
	97	REAL qsurf(ngridmx,nqmx)
	98	REAL q2(ngridmx,nlayermx+1)
	99	! REAL rnaturfi(ngridmx)
	100	real alb(iip1,jjp1),albfi(ngridmx) ! albedos
	101	real ith(iip1,jjp1,nsoilmx),ithfi(ngridmx,nsoilmx) ! thermal inertia (3D)
	102	real surfith(iip1,jjp1),surfithfi(ngridmx) ! surface thermal inertia (2D)
	103	REAL latfi(ngridmx),lonfi(ngridmx),airefi(ngridmx)
	104
	105	INTEGER i,j,l,isoil,ig,idum
	106	real mugaz ! molar mass of the atmosphere
	107
	108	EXTERNAL iniconst,geopot,inigeom
	109	integer ierr !, nbetat
	110	integer ismin
	111	external ismin
	112
	113	c Variables on the new grid along scalar points
	114	c------------------------------------------------------
	115	! REAL p(iip1,jjp1)
	116	REAL t(iip1,jjp1,llm)
	117	real phisold_newgrid(iip1,jjp1)
	118	REAL :: teta(iip1, jjp1, llm)
	119	REAL :: pk(iip1,jjp1,llm)
	120	REAL :: pkf(iip1,jjp1,llm)
	121	REAL :: ps(iip1, jjp1)
	122	REAL :: masse(iip1,jjp1,llm)
	123	REAL :: xpn,xps,xppn(iim),xpps(iim)
	124	REAL :: p3d(iip1, jjp1, llm+1)
	125	REAL :: beta(iip1,jjp1,llm)
	126	! REAL dteta(ip1jmp1,llm)
	127
	128	c Variable de l'ancienne grille
	129	c------------------------------
	130	real time
	131	real tab_cntrl(100)
	132	real tab_cntrl_bis(100)
	133
	134	c variables diverses
	135	c-------------------
	136	real choix_1
	137	character*80 fichnom
	138	integer Lmodif,iq,thermo
	139	character modif*20
	140	real z_reel(iip1,jjp1)
	141	real tsud,albsud,alb_bb,ith_bb,Tiso
	142	real ptoto,pcap,patm,airetot,ptotn,patmn
	143	! real ssum
	144	character*1 yes
	145	logical :: flagiso=.false. , flagps0=.false.
	146	real val, val2, val3 ! to store temporary variables
	147	real :: iceith=2000 ! thermal inertia of subterranean ice
	148	real :: iceithN,iceithS ! values of thermal inertias in N & S hemispheres
	149	integer iref,jref
	150
	151	INTEGER :: itau
	152
	153	INTEGER :: nq,numvanle
	154	character(len=20) :: txt ! to store some text
	155	integer :: count
	156
	157	! MONS data:
	158	real :: MONS_Hdn(iip1,jjp1) ! Hdn: %WEH=Mass fraction of H2O
	159	real :: MONS_d21(iip1,jjp1) ! ice table "depth" (in kg/m2)
	160	! coefficient to apply to convert d21 to 'true' depth (m)
	161	real :: MONS_coeff
	162	real :: MONS_coeffS ! coeff for southern hemisphere
	163	real :: MONS_coeffN ! coeff for northern hemisphere
	164	! real,parameter :: icedepthmin=1.e-3 ! Ice begins at most at that depth
	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 = 744.499 ! for Mars, instead of 1004.70885 (Earth)
	173	preff = 610. ! for Mars, instead of 101325. (Earth)
	174	pa= 20 ! for Mars, instead of 500 (Earth)
	175
	176	c=======================================================================
	177	c Choice of the start file(s) to use
	178	c=======================================================================
	179
	180	write(,) 'From which kind of files do you want to create new',
	181	. 'start and startfi files'
	182	write(,) ' 0 - from a file start_archive'
	183	write(,) ' 1 - from files start and startfi'
	184
	185	c-----------------------------------------------------------------------
	186	c Open file(s) to modify (start or start_archive)
	187	c-----------------------------------------------------------------------
	188
	189	DO
	190	read(,,iostat=ierr) choix_1
	191	if ((choix_1 /= 0).OR.(choix_1 /=1)) EXIT
	192	ENDDO
	193
	194	c Open start_archive
	195	c ~~~~~~~~~~~~~~~~~~~~~~~~~~
	196	if (choix_1.eq.0) then
	197
	198	write(,) 'Creating start files from:'
	199	write(,) './start_archive.nc'
	200	write(,)
	201	fichnom = 'start_archive.nc'
	202	ierr = NF_OPEN (fichnom, NF_NOWRITE,nid)
	203	IF (ierr.NE.NF_NOERR) THEN
	204	write(6,*)' Problem opening file:',fichnom
	205	write(6,*)' ierr = ', ierr
	206	CALL ABORT
	207	ENDIF
	208	tab0 = 50
	209	Lmodif = 1
	210
	211	c OR open start and startfi files
	212	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	213	else
	214	write(,) 'Creating start files from:'
	215	write(,) './start.nc and ./startfi.nc'
	216	write(,)
	217	fichnom = 'start.nc'
	218	ierr = NF_OPEN (fichnom, NF_NOWRITE,nid_dyn)
	219	IF (ierr.NE.NF_NOERR) THEN
	220	write(6,*)' Problem opening file:',fichnom
	221	write(6,*)' ierr = ', ierr
	222	CALL ABORT
	223	ENDIF
	224
	225	fichnom = 'startfi.nc'
	226	ierr = NF_OPEN (fichnom, NF_NOWRITE,nid_fi)
	227	IF (ierr.NE.NF_NOERR) THEN
	228	write(6,*)' Problem opening file:',fichnom
	229	write(6,*)' ierr = ', ierr
	230	CALL ABORT
	231	ENDIF
	232
	233	tab0 = 0
	234	Lmodif = 0
	235
	236	endif
	237
	238	c-----------------------------------------------------------------------
	239	c Lecture du tableau des parametres du run (pour la dynamique)
	240	c-----------------------------------------------------------------------
	241
	242	if (choix_1.eq.0) then
	243
	244	write(,) 'reading tab_cntrl START_ARCHIVE'
	245	c
	246	ierr = NF_INQ_VARID (nid, "controle", nvarid)
	247	#ifdef NC_DOUBLE
	248	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tab_cntrl)
	249	#else
	250	ierr = NF_GET_VAR_REAL(nid, nvarid, tab_cntrl)
	251	#endif
	252	c
	253	else if (choix_1.eq.1) then
	254
	255	write(,) 'reading tab_cntrl START'
	256	c
	257	ierr = NF_INQ_VARID (nid_dyn, "controle", nvarid)
	258	#ifdef NC_DOUBLE
	259	ierr = NF_GET_VAR_DOUBLE(nid_dyn, nvarid, tab_cntrl)
	260	#else
	261	ierr = NF_GET_VAR_REAL(nid_dyn, nvarid, tab_cntrl)
	262	#endif
	263	c
	264	write(,) 'reading tab_cntrl STARTFI'
	265	c
	266	ierr = NF_INQ_VARID (nid_fi, "controle", nvarid)
	267	#ifdef NC_DOUBLE
	268	ierr = NF_GET_VAR_DOUBLE(nid_fi, nvarid, tab_cntrl_bis)
	269	#else
	270	ierr = NF_GET_VAR_REAL(nid_fi, nvarid, tab_cntrl_bis)
	271	#endif
	272	c
	273	do i=1,50
	274	tab_cntrl(i+50)=tab_cntrl_bis(i)
	275	enddo
	276	write(,) 'printing tab_cntrl', tab_cntrl
	277	do i=1,100
	278	write(,) i,tab_cntrl(i)
	279	enddo
	280
	281	endif
	282	c-----------------------------------------------------------------------
	283	c Initialisation des constantes dynamique
	284	c-----------------------------------------------------------------------
	285
	286	kappa = tab_cntrl(9)
	287	etot0 = tab_cntrl(12)
	288	ptot0 = tab_cntrl(13)
	289	ztot0 = tab_cntrl(14)
	290	stot0 = tab_cntrl(15)
	291	ang0 = tab_cntrl(16)
	292	write(,) "Newstart: kappa,etot0,ptot0,ztot0,stot0,ang0"
	293	write(,) kappa,etot0,ptot0,ztot0,stot0,ang0
	294
	295	c-----------------------------------------------------------------------
	296	c Lecture du tab_cntrl et initialisation des constantes physiques
	297	c - pour start: Lmodif = 0 => pas de modifications possibles
	298	c (modif dans le tabfi de readfi + loin)
	299	c - pour start_archive: Lmodif = 1 => modifications possibles
	300	c-----------------------------------------------------------------------
	301	if (choix_1.eq.0) then
	302	call tabfi (nid,Lmodif,tab0,day_ini,lllm,p_rad,
	303	. p_omeg,p_g,p_mugaz,p_daysec,time)
	304	else if (choix_1.eq.1) then
	305	call tabfi (nid_fi,Lmodif,tab0,day_ini,lllm,p_rad,
	306	. p_omeg,p_g,p_mugaz,p_daysec,time)
	307	endif
	308
	309	rad = p_rad
	310	omeg = p_omeg
	311	g = p_g
	312	mugaz = p_mugaz
	313	daysec = p_daysec
	314	! write(,) 'aire',aire
	315
	316
	317	c=======================================================================
	318	c INITIALISATIONS DIVERSES
	319	c=======================================================================
	320	! Load tracer names:
	321	call iniadvtrac(nq,numvanle)
	322
	323	day_step=180 !?! Note: day_step is a common in "control.h"
	324	CALL defrun_new( 99, .TRUE. )
	325	CALL iniconst
	326	CALL inigeom
	327	idum=-1
	328	idum=0
	329
	330	c Initialisation coordonnees /aires
	331	c -------------------------------
	332	! Note: rlatu(:) and rlonv(:) are commons defined in "comgeom.h"
	333	! rlatu() and rlonv() are given in radians
	334	latfi(1)=rlatu(1)
	335	lonfi(1)=0.
	336	DO j=2,jjm
	337	DO i=1,iim
	338	latfi((j-2)*iim+1+i)=rlatu(j)
	339	lonfi((j-2)*iim+1+i)=rlonv(i)
	340	ENDDO
	341	ENDDO
	342	latfi(ngridmx)=rlatu(jjp1)
	343	lonfi(ngridmx)=0.
	344	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
	345
	346	c=======================================================================
	347	c lecture topographie, albedo, inertie thermique, relief sous-maille
	348	c=======================================================================
	349
	350	if (choix_1.ne.1) then ! pour ne pas avoir besoin du fichier
	351	! surface.dat dans le cas des start
	352
	353	c do while((relief(1:3).ne.'mol').AND.(relief(1:3).ne.'pla'))
	354	c write(,)
	355	c write(,) 'choix du relief (mola,pla)'
	356	c write(,) '(Topographie MGS MOLA, plat)'
	357	c read(*,fmt='(a3)') relief
	358	relief="mola"
	359	c enddo
	360
	361	CALL datareadnc(relief,phis,alb,surfith,zmeaS,zstdS,zsigS,zgamS,
	362	. ztheS)
	363
	364	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
	365	! CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ith,ithfi)
	366	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,surfith,surfithfi)
	367	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
	368	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zmeaS,zmea)
	369	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zstdS,zstd)
	370	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zsigS,zsig)
	371	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zgamS,zgam)
	372	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ztheS,zthe)
	373
	374	endif ! of if (choix_1.ne.1)
	375
	376
	377	c=======================================================================
	378	c Lecture des fichiers (start ou start_archive)
	379	c=======================================================================
	380
	381	if (choix_1.eq.0) then
	382
	383	write(,) 'Reading file START_ARCHIVE'
	384	CALL lect_start_archive(date,tsurf,tsoil,emis,q2,
	385	. t,ucov,vcov,ps,co2ice,teta,phisold_newgrid,q,qsurf,
	386	& surfith,nid)
	387	write(,) "OK, read start_archive file"
	388	! copy soil thermal inertia
	389	ithfi(:,:)=inertiedat(:,:)
	390
	391	ierr= NF_CLOSE(nid)
	392
	393	else if (choix_1.eq.1) then ! c'est l'appel a tabfi de phyeta0 qui
	394	! permet de changer les valeurs du
	395	! tab_cntrl Lmodif=1
	396	tab0=0
	397	Lmodif=1 ! Lmodif set to 1 to allow modifications in phyeta0
	398	write(,) 'Reading file START'
	399	fichnom = 'start.nc'
	400	CALL dynetat0(fichnom,nqmx,vcov,ucov,teta,q,masse,
	401	. ps,phis,time)
	402
	403	write(,) 'Reading file STARTFI'
	404	fichnom = 'startfi.nc'
	405	CALL phyetat0 (fichnom,tab0,Lmodif,nsoilmx,nqmx,
	406	. day_ini,time,
	407	. tsurf,tsoil,emis,q2,qsurf,co2ice)
	408
	409	! copy albedo and soil thermal inertia
	410	do i=1,ngridmx
	411	albfi(i) = albedodat(i)
	412	do j=1,nsoilmx
	413	ithfi(i,j) = inertiedat(i,j)
	414	enddo
	415	! build a surfithfi(:) using 1st layer of ithfi(:), which might
	416	! be neede later on if reinitializing soil thermal inertia
	417	surfithfi(i)=ithfi(i,1)
	418	enddo
	419
	420	else
	421	CALL exit(1)
	422	endif
	423
	424	dtvr = daysec/REAL(day_step)
	425	dtphys = dtvr * REAL(iphysiq)
	426
	427	c=======================================================================
	428	c
	429	c=======================================================================
	430	! If tracer names follow 'old' convention (q01, q02, ...) then
	431	! rename them
	432	count=0
	433	do iq=1,nqmx
	434	txt=" "
	435	write(txt,'(a1,i2.2)') 'q',iq
	436	if (txt.eq.tnom(iq)) then
	437	count=count+1
	438	endif
	439	enddo ! of do iq=1,nqmx
	440
	441	if (count.eq.nqmx) then
	442	write(,) 'Newstart: updating tracer names'
	443	! do things the easy but dirty way:
	444	! 1. call inichim_readcallphys (so that callphys.def is read)
	445	call inichim_readcallphys()
	446	! 2. call initracer to set all new tracer names (in noms(:))
	447	call initracer(qsurf,co2ice)
	448	! 3. copy noms(:) to tnom(:)
	449	tnom(1:nqmx)=noms(1:nqmx)
	450	write(,) 'Newstart: updated tracer names'
	451	endif
	452
	453	c=======================================================================
	454	c
	455	c=======================================================================
	456
	457	do ! infinite loop on list of changes
	458
	459	write(,)
	460	write(,)
	461	write(,) 'List of possible changes :'
	462	write(,) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
	463	write(,)
	464	write(,) 'flat : no topography ("aquaplanet")'
	465	write(,) 'bilball : uniform albedo and thermal inertia'
	466	write(,) 'coldspole : cold subsurface and high albedo at S.pole'
	467	write(,) 'qname : change tracer name'
	468	write(,) 'q=0 : ALL tracer =zero'
	469	write(,) 'q=x : give a specific uniform value to one tracer'
	470	write(,) 'ini_q : tracers initialisation for chemistry, water an
	471	$d ice '
	472	write(,) 'ini_q-H2O : tracers initialisation for chemistry and
	473	$ice '
	474	write(,) 'ini_q-iceH2O : tracers initialisation for chemistry on
	475	$ly '
	476	write(,) 'noglacier : Remove tropical H2O ice if \|lat\|<45'
	477	write(,) 'watercapn : H20 ice on permanent N polar cap '
	478	write(,) 'watercaps : H20 ice on permanent S polar cap '
	479	write(,) 'wetstart : start with a wet atmosphere'
	480	write(,) 'isotherm : Isothermal Temperatures, wind set to zero'
	481	write(,) 'co2ice=0 : remove CO2 polar cap'
	482	write(,) 'ptot : change total pressure'
	483	write(,) 'therm_ini_s: Set soil thermal inertia to reference sur
	484	&face values'
	485	write(,) 'subsoilice_n: Put deep underground ice layer in northe
	486	&rn hemisphere'
	487	write(,) 'subsoilice_s: Put deep underground ice layer in southe
	488	&rn hemisphere'
	489	write(,) 'mons_ice: Put underground ice layer according to MONS-
	490	&derived data'
	491
	492	write(,)
	493	write(,) 'Change to perform ?'
	494	write(,) ' (enter keyword or return to end)'
	495	write(,)
	496
	497	read(*,fmt='(a20)') modif
	498	if (modif(1:1) .eq. ' ') exit ! exit loop on changes
	499
	500	write(,)
	501	write(,) modif(1:lnblnk(modif)) , ' : '
	502
	503	c 'flat : no topography ("aquaplanet")'
	504	c -------------------------------------
	505	if (modif(1:lnblnk(modif)) .eq. 'flat') then
	506	c set topo to zero
	507	CALL initial0(ip1jmp1,z_reel)
	508	CALL multscal(ip1jmp1,z_reel,g,phis)
	509	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
	510	write(,) 'topography set to zero.'
	511	write(,) 'WARNING : the subgrid topography parameters',
	512	& ' were not set to zero ! => set calllott to F'
	513
	514	c Choice for surface pressure
	515	yes=' '
	516	do while ((yes.ne.'y').and.(yes.ne.'n'))
	517	write(,) 'Do you wish to choose homogeneous surface',
	518	& 'pressure (y) or let newstart interpolate ',
	519	& ' the previous field (n)?'
	520	read(*,fmt='(a)') yes
	521	end do
	522	if (yes.eq.'y') then
	523	flagps0=.true.
	524	write(,) 'New value for ps (Pa) ?'
	525	201 read(,,iostat=ierr) patm
	526	if(ierr.ne.0) goto 201
	527	write(,)
	528	write(,) ' new ps everywhere (Pa) = ', patm
	529	write(,)
	530	do j=1,jjp1
	531	do i=1,iip1
	532	ps(i,j)=patm
	533	enddo
	534	enddo
	535	end if
	536
	537	c bilball : albedo, inertie thermique uniforme
	538	c --------------------------------------------
	539	else if (modif(1:lnblnk(modif)) .eq. 'bilball') then
	540	write(,) 'constante albedo and iner.therm:'
	541	write(,) 'New value for albedo (ex: 0.25) ?'
	542	101 read(,,iostat=ierr) alb_bb
	543	if(ierr.ne.0) goto 101
	544	write(,)
	545	write(,) ' uniform albedo (new value):',alb_bb
	546	write(,)
	547
	548	write(,) 'New value for thermal inertia (eg: 247) ?'
	549	102 read(,,iostat=ierr) ith_bb
	550	if(ierr.ne.0) goto 102
	551	write(,) 'uniform thermal inertia (new value):',ith_bb
	552	DO j=1,jjp1
	553	DO i=1,iip1
	554	alb(i,j) = alb_bb ! albedo
	555	do isoil=1,nsoilmx
	556	ith(i,j,isoil) = ith_bb ! thermal inertia
	557	enddo
	558	END DO
	559	END DO
	560	! CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ith,ithfi)
	561	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
	562	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
	563
	564	c coldspole : sous-sol de la calotte sud toujours froid
	565	c -----------------------------------------------------
	566	else if (modif(1:lnblnk(modif)) .eq. 'coldspole') then
	567	write(,)'new value for the subsurface temperature',
	568	& ' beneath the permanent southern polar cap ? (eg: 141 K)'
	569	103 read(,,iostat=ierr) tsud
	570	if(ierr.ne.0) goto 103
	571	write(,)
	572	write(,) ' new value of the subsurface temperature:',tsud
	573	c nouvelle temperature sous la calotte permanente
	574	do l=2,nsoilmx
	575	tsoil(ngridmx,l) = tsud
	576	end do
	577
	578
	579	write(,)'new value for the albedo',
	580	& 'of the permanent southern polar cap ? (eg: 0.75)'
	581	104 read(,,iostat=ierr) albsud
	582	if(ierr.ne.0) goto 104
	583	write(,)
	584
	585	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	586	c Option 1: only the albedo of the pole is modified :
	587	albfi(ngridmx)=albsud
	588	write(,) 'ig=',ngridmx,' albedo perennial cap ',
	589	& albfi(ngridmx)
	590
	591	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	592	c Option 2 A haute resolution : coordonnee de la vrai calotte ~
	593	c DO j=1,jjp1
	594	c DO i=1,iip1
	595	c ig=1+(j-2)*iim +i
	596	c if(j.eq.1) ig=1
	597	c if(j.eq.jjp1) ig=ngridmx
	598	c if ((rlatu(j)*180./pi.lt.-84.).and.
	599	c & (rlatu(j)*180./pi.gt.-91.).and.
	600	c & (rlonv(i)*180./pi.gt.-91.).and.
	601	c & (rlonv(i)*180./pi.lt.0.)) then
	602	cc albedo de la calotte permanente fixe a albsud
	603	c alb(i,j)=albsud
	604	c write(,) 'lat=',rlatu(j)*180./pi,
	605	c & ' lon=',rlonv(i)*180./pi
	606	cc fin de la condition sur les limites de la calotte permanente
	607	c end if
	608	c ENDDO
	609	c ENDDO
	610	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	611
	612	c CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
	613
	614
	615	c ptot : Modification of the total pressure: ice + current atmosphere
	616	c -------------------------------------------------------------------
	617	else if (modif(1:lnblnk(modif)) .eq. 'ptot') then
	618
	619	c calcul de la pression totale glace + atm actuelle
	620	patm=0.
	621	airetot=0.
	622	pcap=0.
	623	DO j=1,jjp1
	624	DO i=1,iim
	625	ig=1+(j-2)*iim +i
	626	if(j.eq.1) ig=1
	627	if(j.eq.jjp1) ig=ngridmx
	628	patm = patm + ps(i,j)*aire(i,j)
	629	airetot= airetot + aire(i,j)
	630	pcap = pcap + aire(i,j)co2ice(ig)g
	631	ENDDO
	632	ENDDO
	633	ptoto = pcap + patm
	634
	635	print*,'Current total pressure at surface (co2 ice + atm) ',
	636	& ptoto/airetot
	637
	638	print*,'new value?'
	639	read(,) ptotn
	640	ptotn=ptotn*airetot
	641	patmn=ptotn-pcap
	642	print*,'ptoto,patm,ptotn,patmn'
	643	print*,ptoto,patm,ptotn,patmn
	644	print*,'Mult. factor for pressure (atm only)', patmn/patm
	645	do j=1,jjp1
	646	do i=1,iip1
	647	ps(i,j)=ps(i,j)*patmn/patm
	648	enddo
	649	enddo
	650
	651	c Correction pour la conservation des traceurs
	652	yes=' '
	653	do while ((yes.ne.'y').and.(yes.ne.'n'))
	654	write(,) 'Do you wish to conserve tracer total mass (y)',
	655	& ' or tracer mixing ratio (n) ?'
	656	read(*,fmt='(a)') yes
	657	end do
	658
	659	if (yes.eq.'y') then
	660	write(,) 'OK : conservation of tracer total mass'
	661	DO iq =1, nqmx
	662	DO l=1,llm
	663	DO j=1,jjp1
	664	DO i=1,iip1
	665	q(i,j,l,iq)=q(i,j,l,iq)*patm/patmn
	666	ENDDO
	667	ENDDO
	668	ENDDO
	669	ENDDO
	670	else
	671	write(,) 'OK : conservation of tracer mixing ratio'
	672	end if
	673
	674	c qname : change tracer name
	675	c --------------------------
	676	else if (trim(modif).eq.'qname') then
	677	yes='y'
	678	do while (yes.eq.'y')
	679	write(,) 'Which tracer name do you want to change ?'
	680	do iq=1,nqmx
	681	write(*,'(i3,a3,a20)')iq,' : ',trim(tnom(iq))
	682	enddo
	683	write(*,'(a35,i3)')
	684	& '(enter tracer number; between 1 and ',nqmx
	685	write(,)' or any other value to quit this option)'
	686	read(,) iq
	687	if ((iq.ge.1).and.(iq.le.nqmx)) then
	688	write(,)'Change tracer name ',trim(tnom(iq)),' to ?'
	689	read(,) txt
	690	tnom(iq)=txt
	691	write(,)'Do you want to change another tracer name (y/n)?'
	692	read(*,'(a)') yes
	693	else
	694	! inapropiate value of iq; quit this option
	695	yes='n'
	696	endif ! of if ((iq.ge.1).and.(iq.le.nqmx))
	697	enddo ! of do while (yes.ne.'y')
	698
	699	c q=0 : set tracers to zero
	700	c -------------------------
	701	else if (modif(1:lnblnk(modif)) .eq. 'q=0') then
	702	c mise a 0 des q (traceurs)
	703	write(,) 'Tracers set to 0 (1.E-30 in fact)'
	704	DO iq =1, nqmx
	705	DO l=1,llm
	706	DO j=1,jjp1
	707	DO i=1,iip1
	708	q(i,j,l,iq)=1.e-30
	709	ENDDO
	710	ENDDO
	711	ENDDO
	712	ENDDO
	713
	714	c set surface tracers to zero
	715	DO iq =1, nqmx
	716	DO ig=1,ngridmx
	717	qsurf(ig,iq)=0.
	718	ENDDO
	719	ENDDO
	720
	721	c q=x : initialise tracer manually
	722	c --------------------------------
	723	else if (modif(1:lnblnk(modif)) .eq. 'q=x') then
	724	write(,) 'Which tracer do you want to modify ?'
	725	do iq=1,nqmx
	726	write(,)iq,' : ',trim(tnom(iq))
	727	enddo
	728	write(,) '(choose between 1 and ',nqmx,')'
	729	read(,) iq
	730	write(,)'mixing ratio of tracer ',trim(tnom(iq)),
	731	& ' ? (kg/kg)'
	732	read(,) val
	733	DO l=1,llm
	734	DO j=1,jjp1
	735	DO i=1,iip1
	736	q(i,j,l,iq)=val
	737	ENDDO
	738	ENDDO
	739	ENDDO
	740	write(,) 'SURFACE value of tracer ',trim(tnom(iq)),
	741	& ' ? (kg/m2)'
	742	read(,) val
	743	DO ig=1,ngridmx
	744	qsurf(ig,iq)=val
	745	ENDDO
	746
	747	c ini_q : Initialize tracers for chemistry
	748	c -----------------------------------------------
	749	else if (modif(1:lnblnk(modif)) .eq. 'ini_q') then
	750	c For more than 32 layers, possible to initiate thermosphere only
	751	thermo=0
	752	yes=' '
	753	if (llm.gt.32) then
	754	do while ((yes.ne.'y').and.(yes.ne.'n'))
	755	write(,)'',
	756	& 'initialisation for thermosphere only? (y/n)'
	757	read(*,fmt='(a)') yes
	758	if (yes.eq.'y') then
	759	thermo=1
	760	else
	761	thermo=0
	762	endif
	763	enddo
	764	endif
	765
	766	call inichim_newstart(q,ps,sig,nqmx,latfi,lonfi,airefi,
	767	$ thermo,qsurf)
	768	write(,) 'Chemical species initialized'
	769
	770	if (thermo.eq.0) then
	771	c mise a 0 des qsurf (traceurs a la surface)
	772	DO iq =1, nqmx
	773	DO ig=1,ngridmx
	774	qsurf(ig,iq)=0.
	775	ENDDO
	776	ENDDO
	777	endif
	778
	779	c ini_q-H2O : as above exept for the water vapour tracer
	780	c ------------------------------------------------------
	781	else if (modif(1:lnblnk(modif)) .eq. 'ini_q-H2O') then
	782	! for more than 32 layers, possible to initiate thermosphere only
	783	thermo=0
	784	yes=' '
	785	if(llm.gt.32) then
	786	do while ((yes.ne.'y').and.(yes.ne.'n'))
	787	write(,)'',
	788	& 'initialisation for thermosphere only? (y/n)'
	789	read(*,fmt='(a)') yes
	790	if (yes.eq.'y') then
	791	thermo=1
	792	else
	793	thermo=0
	794	endif
	795	enddo
	796	endif
	797	call inichim_newstart(q,ps,sig,nqmx-1,latfi,lonfi,airefi,
	798	$ thermo,qsurf)
	799	write(,) 'Initialized chem. species exept last (H2O)'
	800
	801	if (thermo.eq.0) then
	802	c set surface tracers to zero, except water ice
	803	DO iq =1, nqmx
	804	if (iq.ne.igcm_h2o_ice) then
	805	DO ig=1,ngridmx
	806	qsurf(ig,iq)=0.
	807	ENDDO
	808	endif
	809	ENDDO
	810	endif
	811
	812	c ini_q-iceH2O : as above exept for ice et H2O
	813	c -----------------------------------------------
	814	else if (modif(1:lnblnk(modif)) .eq. 'ini_q-iceH2O') then
	815	c For more than 32 layers, possible to initiate thermosphere only
	816	thermo=0
	817	yes=' '
	818	if(llm.gt.32) then
	819	do while ((yes.ne.'y').and.(yes.ne.'n'))
	820	write(,)'',
	821	& 'initialisation for thermosphere only? (y/n)'
	822	read(*,fmt='(a)') yes
	823	if (yes.eq.'y') then
	824	thermo=1
	825	else
	826	thermo=0
	827	endif
	828	enddo
	829	endif
	830
	831	call inichim_newstart(q,ps,sig,nqmx-2,latfi,lonfi,airefi,
	832	$ thermo,qsurf)
	833	write(,) 'Initialized chem. species exept ice and H2O'
	834
	835	if (thermo.eq.0) then
	836	c set surface tracers to zero, except water ice
	837	DO iq =1, nqmx
	838	if (iq.ne.igcm_h2o_ice) then
	839	DO ig=1,ngridmx
	840	qsurf(ig,iq)=0.
	841	ENDDO
	842	endif
	843	ENDDO
	844	endif
	845
	846	c wetstart : wet atmosphere with a north to south gradient
	847	c --------------------------------------------------------
	848	else if (modif(1:lnblnk(modif)) .eq. 'wetstart') then
	849	! check that there is indeed a water vapor tracer
	850	if (igcm_h2o_vap.eq.0) then
	851	write(,) "No water vapour tracer! Can't use this option"
	852	stop
	853	endif
	854	DO l=1,llm
	855	DO j=1,jjp1
	856	DO i=1,iip1
	857	q(i,j,l,igcm_h2o_vap)=150.e-6 * (rlatu(j)+pi/2.) / pi
	858	ENDDO
	859	ENDDO
	860	ENDDO
	861
	862	write(,) 'Water mass mixing ratio at north pole='
	863	* ,q(1,1,1,igcm_h2o_vap)
	864	write(,) '---------------------------south pole='
	865	* ,q(1,jjp1,1,igcm_h2o_vap)
	866
	867	c noglacier : remove tropical water ice (to initialize high res sim)
	868	c --------------------------------------------------
	869	else if (modif(1:lnblnk(modif)) .eq. 'noglacier') then
	870	do ig=1,ngridmx
	871	j=(ig-2)/iim +2
	872	if(ig.eq.1) j=1
	873	write(,) 'OK: remove surface ice for \|lat\|<45'
	874	if (abs(rlatu(j)*180./pi).lt.45.) then
	875	qsurf(ig,igcm_h2o_ice)=0.
	876	end if
	877	end do
	878
	879
	880	c watercapn : H20 ice on permanent northern cap
	881	c --------------------------------------------------
	882	else if (modif(1:lnblnk(modif)) .eq. 'watercapn') then
	883	do ig=1,ngridmx
	884	j=(ig-2)/iim +2
	885	if(ig.eq.1) j=1
	886	if (rlatu(j)*180./pi.gt.80.) then
	887	qsurf(ig,igcm_h2o_ice)=1.e5
	888	write(,) 'ig=',ig,' H2O ice mass (kg/m2)= ',
	889	& qsurf(ig,igcm_h2o_ice)
	890	write(,)' ==> Ice mesh South boundary (deg)= ',
	891	& rlatv(j)*180./pi
	892	end if
	893	enddo
	894
	895	c watercaps : H20 ice on permanent southern cap
	896	c -------------------------------------------------
	897	else if (modif(1:lnblnk(modif)) .eq. 'watercaps') then
	898	do ig=1,ngridmx
	899	j=(ig-2)/iim +2
	900	if(ig.eq.1) j=1
	901	if (rlatu(j)*180./pi.lt.-80.) then
	902	qsurf(ig,igcm_h2o_ice)=1.e5
	903	write(,) 'ig=',ig,' H2O ice mass (kg/m2)= ',
	904	& qsurf(ig,igcm_h2o_ice)
	905	write(,)' ==> Ice mesh North boundary (deg)= ',
	906	& rlatv(j-1)*180./pi
	907	end if
	908	enddo
	909
	910	c isotherm : Isothermal temperatures and no winds
	911	c ------------------------------------------------
	912	else if (modif(1:lnblnk(modif)) .eq. 'isotherm') then
	913
	914	write(,)'Isothermal temperature of the atmosphere,
	915	& surface and subsurface'
	916	write(,) 'Value of this temperature ? :'
	917	203 read(,,iostat=ierr) Tiso
	918	if(ierr.ne.0) goto 203
	919
	920	do ig=1, ngridmx
	921	tsurf(ig) = Tiso
	922	end do
	923	do l=2,nsoilmx
	924	do ig=1, ngridmx
	925	tsoil(ig,l) = Tiso
	926	end do
	927	end do
	928	flagiso=.true.
	929	call initial0(llm*ip1jmp1,ucov)
	930	call initial0(llm*ip1jm,vcov)
	931	call initial0(ngridmx*(llm+1),q2)
	932
	933	c co2ice=0 : remove CO2 polar ice caps'
	934	c ------------------------------------------------
	935	else if (modif(1:lnblnk(modif)) .eq. 'co2ice=0') then
	936	do ig=1,ngridmx
	937	co2ice(ig)=0
	938	emis(ig)=emis(ngridmx/2)
	939	end do
	940
	941	! therm_ini_s: (re)-set soil thermal inertia to reference surface values
	942	! ----------------------------------------------------------------------
	943
	944	else if (modif(1:lnblnk(modif)).eq.'therm_ini_s') then
	945	! write(,)"surfithfi(1):",surfithfi(1)
	946	do isoil=1,nsoilmx
	947	inertiedat(1:ngridmx,isoil)=surfithfi(1:ngridmx)
	948	enddo
	949	write(,)'OK: Soil thermal inertia has been reset to referenc
	950	&e surface values'
	951	! write(,)"inertiedat(1,1):",inertiedat(1,1)
	952	ithfi(:,:)=inertiedat(:,:)
	953	! recast ithfi() onto ith()
	954	call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
	955	! Check:
	956	! do i=1,iip1
	957	! do j=1,jjp1
	958	! do isoil=1,nsoilmx
	959	! write(77,*) i,j,isoil," ",ith(i,j,isoil)
	960	! enddo
	961	! enddo
	962	! enddo
	963
	964	! subsoilice_n: Put deep ice layer in northern hemisphere soil
	965	! ------------------------------------------------------------
	966
	967	else if (modif(1:lnblnk(modif)).eq.'subsoilice_n') then
	968
	969	write(,)'From which latitude (in deg.), up to the north pole,
	970	&should we put subterranean ice?'
	971	ierr=1
	972	do while (ierr.ne.0)
	973	read(,,iostat=ierr) val
	974	if (ierr.eq.0) then ! got a value
	975	! do a sanity check
	976	if((val.lt.0.).or.(val.gt.90)) then
	977	write(,)'Latitude should be between 0 and 90 deg. !!!'
	978	ierr=1
	979	else ! find corresponding jref (nearest latitude)
	980	! note: rlatu(:) contains decreasing values of latitude
	981	! starting from PI/2 to -PI/2
	982	do j=1,jjp1
	983	if ((rlatu(j)*180./pi.ge.val).and.
	984	& (rlatu(j+1)*180./pi.le.val)) then
	985	! find which grid point is nearest to val:
	986	if (abs(rlatu(j)*180./pi-val).le.
	987	& abs((rlatu(j+1)*180./pi-val))) then
	988	jref=j
	989	else
	990	jref=j+1
	991	endif
	992
	993	write(,)'Will use nearest grid latitude which is:',
	994	& rlatu(jref)*180./pi
	995	endif
	996	enddo ! of do j=1,jjp1
	997	endif ! of if((val.lt.0.).or.(val.gt.90))
	998	endif !of if (ierr.eq.0)
	999	enddo ! of do while
	1000
	1001	! Build layers() (as in soil_settings.F)
	1002	val2=sqrt(mlayer(0)*mlayer(1))
	1003	val3=mlayer(1)/mlayer(0)
	1004	do isoil=1,nsoilmx
	1005	layer(isoil)=val2(val3*(isoil-1))
	1006	enddo
	1007
	1008	write(,)'At which depth (in m.) does the ice layer begin?'
	1009	write(,)'(currently, the deepest soil layer extends down to:'
	1010	& ,layer(nsoilmx),')'
	1011	ierr=1
	1012	do while (ierr.ne.0)
	1013	read(,,iostat=ierr) val2
	1014	! write(,)'val2:',val2,'ierr=',ierr
	1015	if (ierr.eq.0) then ! got a value, but do a sanity check
	1016	if(val2.gt.layer(nsoilmx)) then
	1017	write(,)'Depth should be less than ',layer(nsoilmx)
	1018	ierr=1
	1019	endif
	1020	if(val2.lt.layer(1)) then
	1021	write(,)'Depth should be more than ',layer(1)
	1022	ierr=1
	1023	endif
	1024	endif
	1025	enddo ! of do while
	1026
	1027	! find the reference index iref the depth corresponds to
	1028	! if (val2.lt.layer(1)) then
	1029	! iref=1
	1030	! else
	1031	do isoil=1,nsoilmx-1
	1032	if((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1)))
	1033	& then
	1034	iref=isoil
	1035	exit
	1036	endif
	1037	enddo
	1038	! endif
	1039
	1040	! write(,)'iref:',iref,' jref:',jref
	1041	! write(,)'layer',layer
	1042	! write(,)'mlayer',mlayer
	1043
	1044	! thermal inertia of the ice:
	1045	ierr=1
	1046	do while (ierr.ne.0)
	1047	write(,)'What is the value of subterranean ice thermal inert
	1048	&ia? (e.g.: 2000)'
	1049	read(,,iostat=ierr)iceith
	1050	enddo ! of do while
	1051
	1052	! recast ithfi() onto ith()
	1053	call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
	1054
	1055	do j=1,jref
	1056	! write(,)'j:',j,'rlatu(j)180./pi:',rlatu(j)180./pi
	1057	do i=1,iip1 ! loop on longitudes
	1058	! Build "equivalent" thermal inertia for the mixed layer
	1059	ith(i,j,iref+1)=sqrt((layer(iref+1)-layer(iref))/
	1060	& (((val2-layer(iref))/(ith(i,j,iref)**2))+
	1061	& ((layer(iref+1)-val2)/(iceith)**2)))
	1062	! Set thermal inertia of lower layers
	1063	do isoil=iref+2,nsoilmx
	1064	ith(i,j,isoil)=iceith ! ice
	1065	enddo
	1066	enddo ! of do i=1,iip1
	1067	enddo ! of do j=1,jjp1
	1068
	1069
	1070	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
	1071
	1072	! do i=1,nsoilmx
	1073	! write(,)'i:',i,'ithfi(1,i):',ithfi(1,i)
	1074	! enddo
	1075
	1076
	1077	! subsoilice_s: Put deep ice layer in southern hemisphere soil
	1078	! ------------------------------------------------------------
	1079
	1080	else if (modif(1:lnblnk(modif)).eq.'subsoilice_s') then
	1081
	1082	write(,)'From which latitude (in deg.), down to the south pol
	1083	&e, should we put subterranean ice?'
	1084	ierr=1
	1085	do while (ierr.ne.0)
	1086	read(,,iostat=ierr) val
	1087	if (ierr.eq.0) then ! got a value
	1088	! do a sanity check
	1089	if((val.gt.0.).or.(val.lt.-90)) then
	1090	write(,)'Latitude should be between 0 and -90 deg. !!!'
	1091	ierr=1
	1092	else ! find corresponding jref (nearest latitude)
	1093	! note: rlatu(:) contains decreasing values of latitude
	1094	! starting from PI/2 to -PI/2
	1095	do j=1,jjp1
	1096	if ((rlatu(j)*180./pi.ge.val).and.
	1097	& (rlatu(j+1)*180./pi.le.val)) then
	1098	! find which grid point is nearest to val:
	1099	if (abs(rlatu(j)*180./pi-val).le.
	1100	& abs((rlatu(j+1)*180./pi-val))) then
	1101	jref=j
	1102	else
	1103	jref=j+1
	1104	endif
	1105
	1106	write(,)'Will use nearest grid latitude which is:',
	1107	& rlatu(jref)*180./pi
	1108	endif
	1109	enddo ! of do j=1,jjp1
	1110	endif ! of if((val.lt.0.).or.(val.gt.90))
	1111	endif !of if (ierr.eq.0)
	1112	enddo ! of do while
	1113
	1114	! Build layers() (as in soil_settings.F)
	1115	val2=sqrt(mlayer(0)*mlayer(1))
	1116	val3=mlayer(1)/mlayer(0)
	1117	do isoil=1,nsoilmx
	1118	layer(isoil)=val2(val3*(isoil-1))
	1119	enddo
	1120
	1121	write(,)'At which depth (in m.) does the ice layer begin?'
	1122	write(,)'(currently, the deepest soil layer extends down to:'
	1123	& ,layer(nsoilmx),')'
	1124	ierr=1
	1125	do while (ierr.ne.0)
	1126	read(,,iostat=ierr) val2
	1127	! write(,)'val2:',val2,'ierr=',ierr
	1128	if (ierr.eq.0) then ! got a value, but do a sanity check
	1129	if(val2.gt.layer(nsoilmx)) then
	1130	write(,)'Depth should be less than ',layer(nsoilmx)
	1131	ierr=1
	1132	endif
	1133	if(val2.lt.layer(1)) then
	1134	write(,)'Depth should be more than ',layer(1)
	1135	ierr=1
	1136	endif
	1137	endif
	1138	enddo ! of do while
	1139
	1140	! find the reference index iref the depth corresponds to
	1141	do isoil=1,nsoilmx-1
	1142	if((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1)))
	1143	& then
	1144	iref=isoil
	1145	exit
	1146	endif
	1147	enddo
	1148
	1149	! write(,)'iref:',iref,' jref:',jref
	1150
	1151	! thermal inertia of the ice:
	1152	ierr=1
	1153	do while (ierr.ne.0)
	1154	write(,)'What is the value of subterranean ice thermal inert
	1155	&ia? (e.g.: 2000)'
	1156	read(,,iostat=ierr)iceith
	1157	enddo ! of do while
	1158
	1159	! recast ithfi() onto ith()
	1160	call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
	1161
	1162	do j=jref,jjp1
	1163	! write(,)'j:',j,'rlatu(j)180./pi:',rlatu(j)180./pi
	1164	do i=1,iip1 ! loop on longitudes
	1165	! Build "equivalent" thermal inertia for the mixed layer
	1166	ith(i,j,iref+1)=sqrt((layer(iref+1)-layer(iref))/
	1167	& (((val2-layer(iref))/(ith(i,j,iref)**2))+
	1168	& ((layer(iref+1)-val2)/(iceith)**2)))
	1169	! Set thermal inertia of lower layers
	1170	do isoil=iref+2,nsoilmx
	1171	ith(i,j,isoil)=iceith ! ice
	1172	enddo
	1173	enddo ! of do i=1,iip1
	1174	enddo ! of do j=jref,jjp1
	1175
	1176
	1177	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
	1178
	1179	c 'mons_ice' : use MONS data to build subsurface ice table
	1180	c --------------------------------------------------------
	1181	else if (modif(1:lnblnk(modif)).eq.'mons_ice') then
	1182
	1183	! 1. Load MONS data
	1184	call load_MONS_data(MONS_Hdn,MONS_d21)
	1185
	1186	! 2. Get parameters from user
	1187	ierr=1
	1188	do while (ierr.ne.0)
	1189	write(,) "Coefficient to apply to MONS 'depth' in Northern",
	1190	& " Hemisphere?"
	1191	write(,) " (should be somewhere between 3.2e-4 and 1.3e-3)"
	1192	read(,,iostat=ierr) MONS_coeffN
	1193	enddo
	1194	ierr=1
	1195	do while (ierr.ne.0)
	1196	write(,) "Coefficient to apply to MONS 'depth' in Southern",
	1197	& " Hemisphere?"
	1198	write(,) " (should be somewhere between 3.2e-4 and 1.3e-3)"
	1199	read(,,iostat=ierr) MONS_coeffS
	1200	enddo
	1201	ierr=1
	1202	do while (ierr.ne.0)
	1203	write(,) "Value of subterranean ice thermal inertia ",
	1204	& " in Northern hemisphere?"
	1205	write(,) " (e.g.: 2000, or perhaps 2290)"
	1206	! read(,,iostat=ierr) iceith
	1207	read(,,iostat=ierr) iceithN
	1208	enddo
	1209	ierr=1
	1210	do while (ierr.ne.0)
	1211	write(,) "Value of subterranean ice thermal inertia ",
	1212	& " in Southern hemisphere?"
	1213	write(,) " (e.g.: 2000, or perhaps 2290)"
	1214	! read(,,iostat=ierr) iceith
	1215	read(,,iostat=ierr) iceithS
	1216	enddo
	1217
	1218	! 3. Build subterranean thermal inertia
	1219
	1220	! initialise subsurface inertia with reference surface values
	1221	do isoil=1,nsoilmx
	1222	ithfi(1:ngridmx,isoil)=surfithfi(1:ngridmx)
	1223	enddo
	1224	! recast ithfi() onto ith()
	1225	call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
	1226
	1227	do i=1,iip1 ! loop on longitudes
	1228	do j=1,jjp1 ! loop on latitudes
	1229	! set MONS_coeff
	1230	if (rlatu(j).ge.0) then ! northern hemisphere
	1231	! N.B: rlatu(:) contains decreasing values of latitude
	1232	! starting from PI/2 to -PI/2
	1233	MONS_coeff=MONS_coeffN
	1234	iceith=iceithN
	1235	else ! southern hemisphere
	1236	MONS_coeff=MONS_coeffS
	1237	iceith=iceithS
	1238	endif
	1239	! check if we should put subterranean ice
	1240	if (MONS_Hdn(i,j).ge.14.0) then ! no ice if Hdn<14%
	1241	! compute depth at which ice lies:
	1242	val=MONS_d21(i,j)*MONS_coeff
	1243	! compute val2= the diurnal skin depth of surface inertia
	1244	! assuming a volumetric heat cap. of C=1.e6 J.m-3.K-1
	1245	val2=ith(i,j,1)1.e-6sqrt(88775./3.14159)
	1246	if (val.lt.val2) then
	1247	! ice must be below the (surface inertia) diurnal skin depth
	1248	val=val2
	1249	endif
	1250	if (val.lt.layer(nsoilmx)) then ! subterranean ice
	1251	! find the reference index iref that depth corresponds to
	1252	iref=0
	1253	do isoil=1,nsoilmx-1
	1254	if ((val.ge.layer(isoil)).and.(val.lt.layer(isoil+1)))
	1255	& then
	1256	iref=isoil
	1257	exit
	1258	endif
	1259	enddo
	1260	! Build "equivalent" thermal inertia for the mixed layer
	1261	ith(i,j,iref+1)=sqrt((layer(iref+1)-layer(iref))/
	1262	& (((val-layer(iref))/(ith(i,j,iref+1)**2))+
	1263	& ((layer(iref+1)-val)/(iceith)**2)))
	1264	! Set thermal inertia of lower layers
	1265	do isoil=iref+2,nsoilmx
	1266	ith(i,j,isoil)=iceith
	1267	enddo
	1268	endif ! of if (val.lt.layer(nsoilmx))
	1269	endif ! of if (MONS_Hdn(i,j).lt.14.0)
	1270	enddo ! do j=1,jjp1
	1271	enddo ! do i=1,iip1
	1272
	1273	! Check:
	1274	! do i=1,iip1
	1275	! do j=1,jjp1
	1276	! do isoil=1,nsoilmx
	1277	! write(77,*) i,j,isoil," ",ith(i,j,isoil)
	1278	! enddo
	1279	! enddo
	1280	! enddo
	1281
	1282	! recast ith() into ithfi()
	1283	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
	1284
	1285	else
	1286	write(,) ' Unknown (misspelled?) option!!!'
	1287	end if ! of if (modif(1:lnblnk(modif)) .eq. '...') elseif ...
	1288
	1289	enddo ! of do ! infinite loop on liste of changes
	1290
	1291	999 continue
	1292
	1293
	1294	c=======================================================================
	1295	c Correct pressure on the new grid (menu 0)
	1296	c=======================================================================
	1297
	1298	if (choix_1.eq.0) then
	1299	r = 1000.*8.31/mugaz
	1300
	1301	do j=1,jjp1
	1302	do i=1,iip1
	1303	ps(i,j) = ps(i,j) *
	1304	. exp((phisold_newgrid(i,j)-phis(i,j)) /
	1305	. (t(i,j,1) * r))
	1306	end do
	1307	end do
	1308
	1309	c periodicity of surface ps in longitude
	1310	do j=1,jjp1
	1311	ps(1,j) = ps(iip1,j)
	1312	end do
	1313	end if
	1314
	1315	c=======================================================================
	1316	c=======================================================================
	1317
	1318	c=======================================================================
	1319	c Initialisation de la physique / ecriture de newstartfi :
	1320	c=======================================================================
	1321
	1322
	1323	CALL inifilr
	1324	CALL pression(ip1jmp1, ap, bp, ps, p3d)
	1325
	1326	c-----------------------------------------------------------------------
	1327	c Initialisation pks:
	1328	c-----------------------------------------------------------------------
	1329
	1330	CALL exner_hyb(ip1jmp1, ps, p3d, beta, pks, pk, pkf)
	1331	! Calcul de la temperature potentielle teta
	1332
	1333	if (flagiso) then
	1334	DO l=1,llm
	1335	DO j=1,jjp1
	1336	DO i=1,iim
	1337	teta(i,j,l) = Tiso * cpp/pk(i,j,l)
	1338	ENDDO
	1339	teta (iip1,j,l)= teta (1,j,l)
	1340	ENDDO
	1341	ENDDO
	1342	else if (choix_1.eq.0) then
	1343	DO l=1,llm
	1344	DO j=1,jjp1
	1345	DO i=1,iim
	1346	teta(i,j,l) = t(i,j,l) * cpp/pk(i,j,l)
	1347	ENDDO
	1348	teta (iip1,j,l)= teta (1,j,l)
	1349	ENDDO
	1350	ENDDO
	1351	endif
	1352
	1353	C Calcul intermediaire
	1354	c
	1355	if (choix_1.eq.0) then
	1356	CALL massdair( p3d, masse )
	1357	c
	1358	print *,' ALPHAX ',alphax
	1359
	1360	DO l = 1, llm
	1361	DO i = 1, iim
	1362	xppn(i) = aire( i, 1 ) * masse( i , 1 , l )
	1363	xpps(i) = aire( i,jjp1 ) * masse( i , jjp1 , l )
	1364	ENDDO
	1365	xpn = SUM(xppn)/apoln
	1366	xps = SUM(xpps)/apols
	1367	DO i = 1, iip1
	1368	masse( i , 1 , l ) = xpn
	1369	masse( i , jjp1 , l ) = xps
	1370	ENDDO
	1371	ENDDO
	1372	endif
	1373	phis(iip1,:) = phis(1,:)
	1374
	1375	CALL inidissip ( lstardis, nitergdiv, nitergrot, niterh,
	1376	* tetagdiv, tetagrot , tetatemp )
	1377	itau=0
	1378	if (choix_1.eq.0) then
	1379	day_ini=int(date)
	1380	endif
	1381	c
	1382	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
	1383
	1384	CALL caldyn0( itau,ucov,vcov,teta,ps,masse,pk,phis ,
	1385	* phi,w, pbaru,pbarv,day_ini+time )
	1386	c CALL caldyn
	1387	c $ ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
	1388	c $ phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, day_ini )
	1389
	1390	CALL dynredem0("restart.nc",day_ini,anneeref,phis,nqmx)
	1391	CALL dynredem1("restart.nc",0.0,vcov,ucov,teta,q,nqmx,masse,ps)
	1392	C
	1393	C Ecriture etat initial physique
	1394	C
	1395
	1396	call physdem1("restartfi.nc",lonfi,latfi,nsoilmx,nqmx,
	1397	. dtphys,real(day_ini),
	1398	. time,tsurf,tsoil,co2ice,emis,q2,qsurf,
	1399	. airefi,albfi,ithfi,zmea,zstd,zsig,zgam,zthe)
	1400
	1401	c=======================================================================
	1402	c Formats
	1403	c=======================================================================
	1404
	1405	1 FORMAT(//10x,'la valeur de im =',i4,2x,'lue sur le fichier de dema
	1406	*rrage est differente de la valeur parametree iim =',i4//)
	1407	2 FORMAT(//10x,'la valeur de jm =',i4,2x,'lue sur le fichier de dema
	1408	*rrage est differente de la valeur parametree jjm =',i4//)
	1409	3 FORMAT(//10x,'la valeur de lllm =',i4,2x,'lue sur le fichier demar
	1410	*rage est differente de la valeur parametree llm =',i4//)
	1411
	1412	end
	1413
	1414	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1415	subroutine load_MONS_data(MONS_Hdn,MONS_d21)
	1416	implicit none
	1417	! routine to load Benedicte Diez MONS dataset, fill in date in southern
	1418	! polar region, and interpolate the result onto the GCM grid
	1419	#include"dimensions.h"
	1420	#include"paramet.h"
	1421	#include"datafile.h"
	1422	#include"comgeom.h"
	1423	! arguments:
	1424	real,intent(out) :: MONS_Hdn(iip1,jjp1) ! Hdn: %WEH=Mass fraction of H2O
	1425	real,intent(out) :: MONS_d21(iip1,jjp1) ! ice table "depth" (in kg/m2)
	1426	! N.B MONS datasets should be of dimension (iip1,jjp1)
	1427	! local variables:
	1428	character(len=88) :: filename="results_MONS_lat_lon_H_depth.txt"
	1429	character(len=88) :: txt ! to store some text
	1430	integer :: ierr,i,j
	1431	integer,parameter :: nblon=180 ! number of longitudes of MONS datasets
	1432	integer,parameter :: nblat=90 ! number of latitudes of MONS datasets
	1433	real :: pi
	1434	real :: longitudes(nblon) ! MONS dataset longitudes
	1435	real :: latitudes(nblat) ! MONS dataset latitudes
	1436	! MONS dataset: mass fraction of H2O where H is assumed to be in H2O
	1437	real :: Hdn(nblon,nblat)
	1438	real :: d21(nblon,nblat)! MONS dataset "depth" (g/cm2)
	1439
	1440	! Extended MONS dataset (for interp_horiz)
	1441	real :: Hdnx(nblon+1,nblat)
	1442	real :: d21x(nblon+1,nblat)
	1443	real :: lon_bound(nblon+1) ! longitude boundaries
	1444	real :: lat_bound(nblat-1) ! latitude boundaries
	1445
	1446	! 1. Initializations:
	1447
	1448	write(,) "Loading MONS data"
	1449
	1450	! Open MONS datafile:
	1451	open(42,file=trim(datafile)//"/"//trim(filename),
	1452	& status="old",iostat=ierr)
	1453	if (ierr/=0) then
	1454	write(,) "Error in load_MONS_data:"
	1455	write(,) "Failed opening file ",
	1456	& trim(datafile)//"/"//trim(filename)
	1457	write(,)'1) You can change the path to the file in '
	1458	write(,)' file phymars/datafile.h'
	1459	write(,)'2) If necessary ',trim(filename),
	1460	& ' (and other datafiles)'
	1461	write(,)' can be obtained online at:'
	1462	write(,)' http://www.lmd.jussieu.fr/~forget/datagcm/datafile'
	1463	CALL ABORT
	1464	else ! skip first line of file (dummy read)
	1465	read(42,*) txt
	1466	endif
	1467
	1468	pi=2.*asin(1.)
	1469
	1470	!2. Load MONS data (on MONS grid)
	1471	do j=1,nblat
	1472	do i=1,nblon
	1473	! swap latitude index so latitudes go from north pole to south pole:
	1474	read(42,*) latitudes(nblat-j+1),longitudes(i),
	1475	& Hdn(i,nblat-j+1),d21(i,nblat-j+1)
	1476	! multiply d21 by 10 to convert from g/cm2 to kg/m2
	1477	d21(i,nblat-j+1)=d21(i,nblat-j+1)*10.0
	1478	enddo
	1479	enddo
	1480	close(42)
	1481
	1482	! there is unfortunately no d21 data for latitudes -77 to -90
	1483	! so we build some by linear interpolation between values at -75
	1484	! and assuming d21=0 at the pole
	1485	do j=84,90 ! latitudes(84)=-77 ; latitudes(83)=-75
	1486	do i=1,nblon
	1487	d21(i,j)=d21(i,83)*((latitudes(j)+90)/15.0)
	1488	enddo
	1489	enddo
	1490
	1491	! 3. Build extended MONS dataset & boundaries (for interp_horiz)
	1492	! longitude boundaries (in radians):
	1493	do i=1,nblon
	1494	! NB: MONS data is every 2 degrees in longitude
	1495	lon_bound(i)=(longitudes(i)+1.0)*pi/180.0
	1496	enddo
	1497	! extra 'modulo' value
	1498	lon_bound(nblon+1)=lon_bound(1)+2.0*pi
	1499
	1500	! latitude boundaries (in radians):
	1501	do j=1,nblat-1
	1502	! NB: Mons data is every 2 degrees in latitude
	1503	lat_bound(j)=(latitudes(j)-1.0)*pi/180.0
	1504	enddo
	1505
	1506	! MONS datasets:
	1507	do j=1,nblat
	1508	Hdnx(1:nblon,j)=Hdn(1:nblon,j)
	1509	Hdnx(nblon+1,j)=Hdnx(1,j)
	1510	d21x(1:nblon,j)=d21(1:nblon,j)
	1511	d21x(nblon+1,j)=d21x(1,j)
	1512	enddo
	1513
	1514	! Interpolate onto GCM grid
	1515	call interp_horiz(Hdnx,MONS_Hdn,nblon,nblat-1,iim,jjm,1,
	1516	& lon_bound,lat_bound,rlonu,rlatv)
	1517	call interp_horiz(d21x,MONS_d21,nblon,nblat-1,iim,jjm,1,
	1518	& lon_bound,lat_bound,rlonu,rlatv)
	1519
	1520	end subroutine

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

Download in other formats:

Original Format