Context Navigation

newstart.F @ 937

Last change on this file since 937 was 697, checked in by emillour, 13 years ago

Mars GCM:

Minor cleanup in surfini.F
Corrected the polar mesh surface area which was wrong in the physics (changes in phyetat0.F, calfis.F and newstart.F)

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

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

Download in other formats:

Original Format