Context Navigation

source: trunk/LMDZ.GENERIC/libf/dyn3d/newstart.old @ 211

Last change on this file since 211 was 135, checked in by aslmd, 14 years ago
CHANGEMENT ARBORESCENCE ETAPE 2 -- NON COMPLET
File size: 51.3 KB

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

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

Download in other formats: