Context Navigation

newstart.F @ 629

Last change on this file since 629 was 618, checked in by emillour, 13 years ago

Mars GCM:

Update of the chemistry package:

calchim.F90 :
- upgraded from calchim.F
- can run with or without CH4
- fixed the mass conservation scheme
photochemistry.F : -chemistry timestep is now independant from physics timestep -can run with or without a CH4 tracer
- removed initial tests on species, since these are already done in calchim
removed inichim_readcallphys.F (not used any more)
concentrations.F: adaptated to better handle indexes and molar masses of

tracers. "ncomp" (in chimiedata.h) no longer needed.

inichim_newstart.F90 : rewritten and cleaned (won't be compatible with

old start files where tracer names are q01,...).
Now handles hybrid levels and automaticaly adapts
depending on which tracers are available.

newstart.F : adapted to followup changes in inchim_newstart.F90 and some

cleanup around the initialization of tracer names and surface
values.

FL+EM

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

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

Download in other formats:

Original Format