Context Navigation

newstart.F @ 480

Last change on this file since 480 was 480, checked in by tnavarro, 13 years ago

new option in newstart: ini_h2osurf, water ice surface initialisation to keep seasonal frost

but get rid of polar cap positive or negative values.

File size: 52.7 KB

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

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

Download in other formats:

Original Format