Context Navigation

newstart.F @ 1422

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

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

Download in other formats:

Original Format