Context Navigation

newstart.F @ 1241

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

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

Download in other formats:

Original Format