Context Navigation

newstart.F @ 607

Last change on this file since 607 was 563, checked in by emillour, 13 years ago

Mars GCM:

Added option "q=profile" in newstart to initialize a given tracer with

a profile specified in an ASCII file.

Minor fix on lect_start_archive : when tracers were not found, a default

value was asked to user, but twice! Once is enough.

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

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

Download in other formats:

Original Format