Context Navigation

newstart.F @ 1047

Last change on this file since 1047 was 1047, checked in by emillour, 11 years ago

Mars GCM:

IMPORTANT CHANGE: Removed all reference/use of ngridmx (dimphys.h) in routines (necessary prerequisite to using parallel dynamics); in most cases this just means adding 'ngrid' as routine argument, and making local saved variables allocatable (and allocated at first call). In the process, had to convert many *.h files to equivalent modules: yomaer.h => yomaer_h.F90 , surfdat.h => surfdat_h.F90 , comsaison.h => comsaison_h.F90 , yomlw.h => yomlw_h.F90 , comdiurn.h => comdiurn_h.F90 , dimradmars.h => dimradmars_mod.F90 , comgeomfi.h => comgeomfi_h.F90, comsoil.h => comsoil_h.F90 , slope.h => slope_mod.F90
Also updated EOF routines, everything is now in eofdump_mod.F90
Removed unused routine lectfux.F (in dyn3d)

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

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

Download in other formats:

Original Format