Context Navigation

newstart.F @ 699

Last change on this file since 699 was 699, checked in by emillour, 13 years ago

Generic GCM:

Corrected the polar mesh surface area which was wrong in the physics (changes in phyetat0.F, calfis.F and newstart.F)
Some cleanup in newstart.F (removed some obsolete "Mars" options: mons_ice,..) and also added option "q=profile" to initialize a tracer with a profile read from file "profile_tracername"

File size: 47.9 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_cpp,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 emisread ! added by RW
	91	REAL qsurf(ngridmx,nqmx)
	92	REAL q2(ngridmx,nlayermx+1)
	93	! REAL rnaturfi(ngridmx)
	94	real alb(iip1,jjp1),albfi(ngridmx) ! albedos
	95	real ith(iip1,jjp1,nsoilmx),ithfi(ngridmx,nsoilmx) ! thermal inertia (3D)
	96	real surfith(iip1,jjp1),surfithfi(ngridmx) ! surface thermal inertia (2D)
	97	REAL latfi(ngridmx),lonfi(ngridmx),airefi(ngridmx)
	98
	99	INTEGER i,j,l,isoil,ig,idum
	100	real mugaz ! molar mass of the atmosphere
	101
	102	integer ierr
	103
	104	c Variables on the new grid along scalar points
	105	c------------------------------------------------------
	106	! REAL p(iip1,jjp1)
	107	REAL t(iip1,jjp1,llm)
	108	REAL tset(iip1,jjp1,llm)
	109	real phisold_newgrid(iip1,jjp1)
	110	REAL :: teta(iip1, jjp1, llm)
	111	REAL :: pk(iip1,jjp1,llm)
	112	REAL :: pkf(iip1,jjp1,llm)
	113	REAL :: ps(iip1, jjp1)
	114	REAL :: masse(iip1,jjp1,llm)
	115	REAL :: xpn,xps,xppn(iim),xpps(iim)
	116	REAL :: p3d(iip1, jjp1, llm+1)
	117	REAL :: beta(iip1,jjp1,llm)
	118	! REAL dteta(ip1jmp1,llm)
	119
	120	c Variable de l'ancienne grille
	121	c------------------------------
	122	real time
	123	real tab_cntrl(100)
	124	real tab_cntrl_bis(100)
	125
	126	c variables diverses
	127	c-------------------
	128	real choix_1,pp
	129	character*80 fichnom
	130	integer Lmodif,iq,thermo
	131	character modif*20
	132	real z_reel(iip1,jjp1)
	133	real tsud,albsud,alb_bb,ith_bb,Tiso,Tabove
	134	real ptoto,pcap,patm,airetot,ptotn,patmn,psea
	135	! real ssum
	136	character*1 yes
	137	logical :: flagtset=.false. , flagps0=.false.
	138	real val, val2, val3 ! to store temporary variables
	139	real :: iceith=2000 ! thermal inertia of subterranean ice
	140	integer iref,jref
	141
	142	INTEGER :: itau
	143
	144	INTEGER :: nq,numvanle
	145	character(len=20) :: txt ! to store some text
	146	integer :: count
[699]	147	real :: profile(llm+1) ! to store an atmospheric profile + surface value
[135]	148
[253]	149	! added by RW for test
	150	real pmean, phi0
	151
	152	! added by BC for equilibrium temperature startup
	153	real teque
	154
	155	! added by BC for cloud fraction setup
	156	REAL hice(ngridmx),cloudfrac(ngridmx,nlayermx)
	157	REAL totalfrac(ngridmx)
	158
	159	! added by RW for nuketharsis
	160	real fact1
	161	real fact2
	162
	163
[135]	164	c sortie visu pour les champs dynamiques
	165	c---------------------------------------
	166	! INTEGER :: visuid
	167	! real :: time_step,t_ops,t_wrt
	168	! CHARACTER*80 :: visu_file
	169
[535]	170	cpp = 0.
	171	preff = 0.
	172	pa = 0. ! to ensure disaster rather than minor error if we don`t
	173	! make deliberate choice of these values elsewhere.
[135]	174
	175	c=======================================================================
	176	c Choice of the start file(s) to use
	177	c=======================================================================
	178
	179	write(,) 'From which kind of files do you want to create new',
	180	. 'start and startfi files'
	181	write(,) ' 0 - from a file start_archive'
	182	write(,) ' 1 - from files start and startfi'
	183
	184	c-----------------------------------------------------------------------
	185	c Open file(s) to modify (start or start_archive)
	186	c-----------------------------------------------------------------------
	187
	188	DO
	189	read(,,iostat=ierr) choix_1
	190	if ((choix_1 /= 0).OR.(choix_1 /=1)) EXIT
	191	ENDDO
	192
	193	c Open start_archive
	194	c ~~~~~~~~~~~~~~~~~~~~~~~~~~
	195	if (choix_1.eq.0) then
	196
	197	write(,) 'Creating start files from:'
	198	write(,) './start_archive.nc'
	199	write(,)
	200	fichnom = 'start_archive.nc'
	201	ierr = NF_OPEN (fichnom, NF_NOWRITE,nid)
	202	IF (ierr.NE.NF_NOERR) THEN
	203	write(6,*)' Problem opening file:',fichnom
	204	write(6,*)' ierr = ', ierr
	205	CALL ABORT
	206	ENDIF
	207	tab0 = 50
	208	Lmodif = 1
	209
	210	c OR open start and startfi files
	211	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	212	else
	213	write(,) 'Creating start files from:'
	214	write(,) './start.nc and ./startfi.nc'
	215	write(,)
	216	fichnom = 'start.nc'
	217	ierr = NF_OPEN (fichnom, NF_NOWRITE,nid_dyn)
	218	IF (ierr.NE.NF_NOERR) THEN
	219	write(6,*)' Problem opening file:',fichnom
	220	write(6,*)' ierr = ', ierr
	221	CALL ABORT
	222	ENDIF
	223
	224	fichnom = 'startfi.nc'
	225	ierr = NF_OPEN (fichnom, NF_NOWRITE,nid_fi)
	226	IF (ierr.NE.NF_NOERR) THEN
	227	write(6,*)' Problem opening file:',fichnom
	228	write(6,*)' ierr = ', ierr
	229	CALL ABORT
	230	ENDIF
	231
	232	tab0 = 0
	233	Lmodif = 0
	234
	235	endif
	236
[649]	237
	238	c=======================================================================
	239	c INITIALISATIONS DIVERSES
	240	c=======================================================================
	241	! Load tracer names:
	242	call iniadvtrac(nq,numvanle)
	243	! tnom(:) now contains tracer names
	244	! JL12 we will need the tracer names to read start in dyneta0
	245
[135]	246	c-----------------------------------------------------------------------
	247	c Lecture du tableau des parametres du run (pour la dynamique)
	248	c-----------------------------------------------------------------------
	249
	250	if (choix_1.eq.0) then
	251
	252	write(,) 'reading tab_cntrl START_ARCHIVE'
	253	c
	254	ierr = NF_INQ_VARID (nid, "controle", nvarid)
	255	#ifdef NC_DOUBLE
	256	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tab_cntrl)
	257	#else
	258	ierr = NF_GET_VAR_REAL(nid, nvarid, tab_cntrl)
	259	#endif
	260	c
	261	else if (choix_1.eq.1) then
	262
	263	write(,) 'reading tab_cntrl START'
	264	c
	265	ierr = NF_INQ_VARID (nid_dyn, "controle", nvarid)
	266	#ifdef NC_DOUBLE
	267	ierr = NF_GET_VAR_DOUBLE(nid_dyn, nvarid, tab_cntrl)
	268	#else
	269	ierr = NF_GET_VAR_REAL(nid_dyn, nvarid, tab_cntrl)
	270	#endif
	271	c
	272	write(,) 'reading tab_cntrl STARTFI'
	273	c
	274	ierr = NF_INQ_VARID (nid_fi, "controle", nvarid)
	275	#ifdef NC_DOUBLE
	276	ierr = NF_GET_VAR_DOUBLE(nid_fi, nvarid, tab_cntrl_bis)
	277	#else
	278	ierr = NF_GET_VAR_REAL(nid_fi, nvarid, tab_cntrl_bis)
	279	#endif
	280	c
	281	do i=1,50
	282	tab_cntrl(i+50)=tab_cntrl_bis(i)
	283	enddo
[649]	284	write(,) 'printing tab_cntrl', tab_cntrl
	285	do i=1,100
	286	write(,) i,tab_cntrl(i)
	287	enddo
	288
	289	! Lmodif set to 0 to disable modifications possibility in phyeta0
	290	write(,) 'Reading file START'
	291	fichnom = 'start.nc'
	292	CALL dynetat0(fichnom,nqmx,vcov,ucov,teta,q,masse,
	293	. ps,phis,time)
	294
	295	write(,) 'Reading file STARTFI'
	296	fichnom = 'startfi.nc'
	297	CALL phyetat0 (fichnom,tab0,Lmodif,nsoilmx,nqmx,
	298	. day_ini,time,
	299	. tsurf,tsoil,emis,q2,qsurf, !) ! temporary modif by RDW
	300	. cloudfrac,totalfrac,hice)
	301
	302	! copy albedo and soil thermal inertia
	303	do i=1,ngridmx
	304	albfi(i) = albedodat(i)
	305	do j=1,nsoilmx
	306	ithfi(i,j) = inertiedat(i,j)
	307	enddo
	308	! build a surfithfi(:) using 1st layer of ithfi(:), which might
	309	! be neede later on if reinitializing soil thermal inertia
	310	surfithfi(i)=ithfi(i,1)
	311	enddo
	312
[135]	313
	314	endif
	315	c-----------------------------------------------------------------------
	316	c Initialisation des constantes dynamique
	317	c-----------------------------------------------------------------------
	318
	319	kappa = tab_cntrl(9)
	320	etot0 = tab_cntrl(12)
	321	ptot0 = tab_cntrl(13)
	322	ztot0 = tab_cntrl(14)
	323	stot0 = tab_cntrl(15)
	324	ang0 = tab_cntrl(16)
	325	write(,) "Newstart: kappa,etot0,ptot0,ztot0,stot0,ang0"
	326	write(,) kappa,etot0,ptot0,ztot0,stot0,ang0
	327
	328	! for vertical coordinate
	329	preff=tab_cntrl(18) ! reference surface pressure
	330	pa=tab_cntrl(17) ! reference pressure at which coord is purely pressure
	331	!NB: in start_archive files tab_cntrl(17)=tab_cntrl(18)=0
	332	write(,) "Newstart: preff=",preff," pa=",pa
	333	yes=' '
	334	do while ((yes.ne.'y').and.(yes.ne.'n'))
	335	write(,) "Change the values of preff and pa? (y/n)"
	336	read(*,fmt='(a)') yes
	337	end do
	338	if (yes.eq.'y') then
	339	write(,)"New value of reference surface pressure preff?"
	340	write(,)" (for Mars, typically preff=610)"
	341	read(,) preff
	342	write(,)"New value of reference pressure pa for purely"
	343	write(,)"pressure levels (for hybrid coordinates)?"
	344	write(,)" (for Mars, typically pa=20)"
	345	read(,) pa
	346	endif
	347	c-----------------------------------------------------------------------
	348	c Lecture du tab_cntrl et initialisation des constantes physiques
	349	c - pour start: Lmodif = 0 => pas de modifications possibles
	350	c (modif dans le tabfi de readfi + loin)
	351	c - pour start_archive: Lmodif = 1 => modifications possibles
	352	c-----------------------------------------------------------------------
	353	if (choix_1.eq.0) then
	354	call tabfi (nid,Lmodif,tab0,day_ini,lllm,p_rad,
	355	. p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
	356	else if (choix_1.eq.1) then
[649]	357	Lmodif=1 ! Lmodif set to 1 to allow modifications in phyeta0
[135]	358	call tabfi (nid_fi,Lmodif,tab0,day_ini,lllm,p_rad,
	359	. p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
	360	endif
	361
	362	rad = p_rad
	363	omeg = p_omeg
	364	g = p_g
	365	cpp = p_cpp
	366	mugaz = p_mugaz
	367	daysec = p_daysec
	368
[535]	369	kappa = 8.3141000./(p_mugazp_cpp) ! added by RDW
[135]	370
	371
	372	c=======================================================================
	373	c INITIALISATIONS DIVERSES
	374	c=======================================================================
	375	! Initialize global tracer indexes (stored in tracer.h)
	376	call initracer()
	377	! Load parameters from run.def file
	378	CALL defrun_new( 99, .TRUE. )
	379	CALL iniconst
	380	CALL inigeom
	381	idum=-1
	382	idum=0
	383
	384	c Initialisation coordonnees /aires
	385	c -------------------------------
	386	! Note: rlatu(:) and rlonv(:) are commons defined in "comgeom.h"
	387	! rlatu() and rlonv() are given in radians
	388	latfi(1)=rlatu(1)
	389	lonfi(1)=0.
	390	DO j=2,jjm
	391	DO i=1,iim
	392	latfi((j-2)*iim+1+i)=rlatu(j)
	393	lonfi((j-2)*iim+1+i)=rlonv(i)
	394	ENDDO
	395	ENDDO
	396	latfi(ngridmx)=rlatu(jjp1)
	397	lonfi(ngridmx)=0.
[699]	398
	399	! build airefi(), mesh area on physics grid
[135]	400	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
[699]	401	! Poles are single points on physics grid
	402	airefi(1)=airefi(1)*iim
	403	airefi(ngridmx)=airefi(ngridmx)*iim
[135]	404
	405	c=======================================================================
	406	c lecture topographie, albedo, inertie thermique, relief sous-maille
	407	c=======================================================================
	408
	409	if (choix_1.ne.1) then ! pour ne pas avoir besoin du fichier
	410	! surface.dat dans le cas des start
	411
	412	c do while((relief(1:3).ne.'mol').AND.(relief(1:3).ne.'pla'))
	413	c write(,)
	414	c write(,) 'choix du relief (mola,pla)'
	415	c write(,) '(Topographie MGS MOLA, plat)'
	416	c read(*,fmt='(a3)') relief
	417	relief="mola"
	418	c enddo
	419
	420	CALL datareadnc(relief,phis,alb,surfith,zmeaS,zstdS,zsigS,zgamS,
	421	. ztheS)
	422
	423	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
	424	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,surfith,surfithfi)
	425	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
	426	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zmeaS,zmea)
	427	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zstdS,zstd)
	428	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zsigS,zsig)
	429	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zgamS,zgam)
	430	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ztheS,zthe)
	431
	432	endif ! of if (choix_1.ne.1)
	433
	434
	435	c=======================================================================
	436	c Lecture des fichiers (start ou start_archive)
	437	c=======================================================================
	438
	439	if (choix_1.eq.0) then
	440
	441	write(,) 'Reading file START_ARCHIVE'
	442	CALL lect_start_archive(date,tsurf,tsoil,emis,q2,
	443	. t,ucov,vcov,ps,teta,phisold_newgrid,q,qsurf,
	444	& surfith,nid)
	445	write(,) "OK, read start_archive file"
	446	! copy soil thermal inertia
	447	ithfi(:,:)=inertiedat(:,:)
	448
	449	ierr= NF_CLOSE(nid)
	450
[649]	451	else if (choix_1.eq.1) then
	452	!do nothing, start and startfi have already been read
[135]	453	else
	454	CALL exit(1)
	455	endif
	456
	457	dtvr = daysec/FLOAT(day_step)
	458	dtphys = dtvr * FLOAT(iphysiq)
	459
	460	c=======================================================================
	461	c
	462	c=======================================================================
	463
	464	do ! infinite loop on list of changes
	465
	466	write(,)
	467	write(,)
	468	write(,) 'List of possible changes :'
	469	write(,) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
	470	write(,)
	471	write(,) 'flat : no topography ("aquaplanet")'
[253]	472	write(,) 'nuketharsis : no Tharsis bulge'
[135]	473	write(,) 'bilball : uniform albedo and thermal inertia'
	474	write(,) 'coldspole : cold subsurface and high albedo at S.pole'
	475	write(,) 'qname : change tracer name'
	476	write(,) 'q=0 : ALL tracer =zero'
	477	write(,) 'q=x : give a specific uniform value to one tracer'
[699]	478	write(,) 'q=profile : specify a profile for a tracer'
	479	! write(,) 'ini_q : tracers initialisation for chemistry, water an
	480	! $d ice '
	481	! write(,) 'ini_q-H2O : tracers initialisation for chemistry and
	482	! $ice '
	483	! write(,) 'ini_q-iceH2O : tracers initialisation for chemistry on
	484	! $ly '
[135]	485	write(,) 'noglacier : Remove tropical H2O ice if \|lat\|<45'
	486	write(,) 'watercapn : H20 ice on permanent N polar cap '
	487	write(,) 'watercaps : H20 ice on permanent S polar cap '
[253]	488	write(,) 'noacglac : H2O ice across Noachis Terra'
[135]	489	write(,) 'oborealis : H2O ice across Vastitas Borealis'
[253]	490	write(,) 'iceball : Thick ice layer all over surface'
[135]	491	write(,) 'wetstart : start with a wet atmosphere'
[253]	492	write(,) 'isotherm : Isothermal Temperatures, wind set to zero'
	493	write(,) 'radequi : Earth-like radiative equilibrium temperature
	494	$ profile (lat-alt) and winds set to zero'
[135]	495	write(,) 'coldstart : Start X K above the CO2 frost point and
	496	$set wind to zero (assumes 100% CO2)'
	497	write(,) 'co2ice=0 : remove CO2 polar cap'
	498	write(,) 'ptot : change total pressure'
	499	write(,) 'emis : change surface emissivity'
[253]	500	write(,) 'therm_ini_s : Set soil thermal inertia to reference sur
[135]	501	&face values'
	502
	503	write(,)
	504	write(,) 'Change to perform ?'
	505	write(,) ' (enter keyword or return to end)'
	506	write(,)
	507
	508	read(*,fmt='(a20)') modif
	509	if (modif(1:1) .eq. ' ') exit ! exit loop on changes
	510
	511	write(,)
	512	write(,) trim(modif) , ' : '
	513
	514	c 'flat : no topography ("aquaplanet")'
	515	c -------------------------------------
[699]	516	if (trim(modif) .eq. 'flat') then
[135]	517	c set topo to zero
[699]	518	z_reel(1:iip1,1:jjp1)=0
	519	phis(1:iip1,1:jjp1)=g*z_reel(1:iip1,1:jjp1)
[135]	520	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
	521	write(,) 'topography set to zero.'
	522	write(,) 'WARNING : the subgrid topography parameters',
	523	& ' were not set to zero ! => set calllott to F'
	524
[253]	525	c Choice of surface pressure
[135]	526	yes=' '
	527	do while ((yes.ne.'y').and.(yes.ne.'n'))
	528	write(,) 'Do you wish to choose homogeneous surface',
	529	& 'pressure (y) or let newstart interpolate ',
	530	& ' the previous field (n)?'
	531	read(*,fmt='(a)') yes
	532	end do
	533	if (yes.eq.'y') then
	534	flagps0=.true.
	535	write(,) 'New value for ps (Pa) ?'
	536	201 read(,,iostat=ierr) patm
	537	if(ierr.ne.0) goto 201
	538	write(,)
	539	write(,) ' new ps everywhere (Pa) = ', patm
	540	write(,)
	541	do j=1,jjp1
	542	do i=1,iip1
	543	ps(i,j)=patm
	544	enddo
	545	enddo
	546	end if
	547
[253]	548	c 'nuketharsis : no tharsis bulge for Early Mars'
	549	c ---------------------------------------------
[699]	550	else if (trim(modif) .eq. 'nuketharsis') then
[253]	551
	552	DO j=1,jjp1
	553	DO i=1,iim
	554	ig=1+(j-2)*iim +i
	555	if(j.eq.1) ig=1
	556	if(j.eq.jjp1) ig=ngridmx
	557
	558	fact1=(((rlonv(i)180./pi)+100)*2 +
	559	& (rlatu(j)180./pi)2)/65*2
	560	fact2=exp( -fact1**2.5 )
	561
	562	phis(i,j) = phis(i,j) - (phis(i,j)+4000.g)fact2
	563
	564	! if(phis(i,j).gt.2500.*g)then
	565	! if(rlatu(j)*180./pi.gt.-80.)then ! avoid chopping south polar cap
	566	! phis(i,j)=2500.*g
	567	! endif
	568	! endif
	569
	570	ENDDO
	571	ENDDO
	572	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
	573
	574
	575	c bilball : uniform albedo, thermal inertia
	576	c -----------------------------------------
[699]	577	else if (trim(modif) .eq. 'bilball') then
[135]	578	write(,) 'constante albedo and iner.therm:'
	579	write(,) 'New value for albedo (ex: 0.25) ?'
	580	101 read(,,iostat=ierr) alb_bb
	581	if(ierr.ne.0) goto 101
	582	write(,)
	583	write(,) ' uniform albedo (new value):',alb_bb
	584	write(,)
	585
	586	write(,) 'New value for thermal inertia (eg: 247) ?'
	587	102 read(,,iostat=ierr) ith_bb
	588	if(ierr.ne.0) goto 102
	589	write(,) 'uniform thermal inertia (new value):',ith_bb
	590	DO j=1,jjp1
	591	DO i=1,iip1
	592	alb(i,j) = alb_bb ! albedo
	593	do isoil=1,nsoilmx
	594	ith(i,j,isoil) = ith_bb ! thermal inertia
	595	enddo
	596	END DO
	597	END DO
	598	! CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ith,ithfi)
	599	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
	600	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
	601
	602	c coldspole : sous-sol de la calotte sud toujours froid
	603	c -----------------------------------------------------
[699]	604	else if (trim(modif) .eq. 'coldspole') then
[135]	605	write(,)'new value for the subsurface temperature',
	606	& ' beneath the permanent southern polar cap ? (eg: 141 K)'
	607	103 read(,,iostat=ierr) tsud
	608	if(ierr.ne.0) goto 103
	609	write(,)
	610	write(,) ' new value of the subsurface temperature:',tsud
	611	c nouvelle temperature sous la calotte permanente
	612	do l=2,nsoilmx
	613	tsoil(ngridmx,l) = tsud
	614	end do
	615
	616
	617	write(,)'new value for the albedo',
	618	& 'of the permanent southern polar cap ? (eg: 0.75)'
	619	104 read(,,iostat=ierr) albsud
	620	if(ierr.ne.0) goto 104
	621	write(,)
	622
	623	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	624	c Option 1: only the albedo of the pole is modified :
	625	albfi(ngridmx)=albsud
	626	write(,) 'ig=',ngridmx,' albedo perennial cap ',
	627	& albfi(ngridmx)
	628
	629	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	630	c Option 2 A haute resolution : coordonnee de la vrai calotte ~
	631	c DO j=1,jjp1
	632	c DO i=1,iip1
	633	c ig=1+(j-2)*iim +i
	634	c if(j.eq.1) ig=1
	635	c if(j.eq.jjp1) ig=ngridmx
	636	c if ((rlatu(j)*180./pi.lt.-84.).and.
	637	c & (rlatu(j)*180./pi.gt.-91.).and.
	638	c & (rlonv(i)*180./pi.gt.-91.).and.
	639	c & (rlonv(i)*180./pi.lt.0.)) then
	640	cc albedo de la calotte permanente fixe a albsud
	641	c alb(i,j)=albsud
	642	c write(,) 'lat=',rlatu(j)*180./pi,
	643	c & ' lon=',rlonv(i)*180./pi
	644	cc fin de la condition sur les limites de la calotte permanente
	645	c end if
	646	c ENDDO
	647	c ENDDO
	648	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	649
	650	c CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
	651
	652
	653	c ptot : Modification of the total pressure: ice + current atmosphere
	654	c -------------------------------------------------------------------
[699]	655	else if (trim(modif).eq.'ptot') then
[135]	656
	657	! check if we have a co2_ice surface tracer:
	658	if (igcm_co2_ice.eq.0) then
	659	write(,) " No surface CO2 ice !"
	660	write(,) " only atmospheric pressure will be considered!"
	661	endif
	662	c calcul de la pression totale glace + atm actuelle
	663	patm=0.
	664	airetot=0.
	665	pcap=0.
	666	DO j=1,jjp1
	667	DO i=1,iim
	668	ig=1+(j-2)*iim +i
	669	if(j.eq.1) ig=1
	670	if(j.eq.jjp1) ig=ngridmx
	671	patm = patm + ps(i,j)*aire(i,j)
	672	airetot= airetot + aire(i,j)
	673	if (igcm_co2_ice.ne.0) then
	674	!pcap = pcap + aire(i,j)co2ice(ig)g
	675	pcap = pcap + aire(i,j)qsurf(ig,igcm_co2_ice)g
	676	endif
	677	ENDDO
	678	ENDDO
	679	ptoto = pcap + patm
	680
	681	print*,'Current total pressure at surface (co2 ice + atm) ',
	682	& ptoto/airetot
	683
	684	print*,'new value?'
	685	read(,) ptotn
	686	ptotn=ptotn*airetot
	687	patmn=ptotn-pcap
	688	print*,'ptoto,patm,ptotn,patmn'
	689	print*,ptoto,patm,ptotn,patmn
	690	print*,'Mult. factor for pressure (atm only)', patmn/patm
	691	do j=1,jjp1
	692	do i=1,iip1
	693	ps(i,j)=ps(i,j)*patmn/patm
	694	enddo
	695	enddo
	696
[253]	697
	698
[135]	699	c Correction pour la conservation des traceurs
	700	yes=' '
	701	do while ((yes.ne.'y').and.(yes.ne.'n'))
	702	write(,) 'Do you wish to conserve tracer total mass (y)',
	703	& ' or tracer mixing ratio (n) ?'
	704	read(*,fmt='(a)') yes
	705	end do
	706
	707	if (yes.eq.'y') then
	708	write(,) 'OK : conservation of tracer total mass'
	709	DO iq =1, nqmx
	710	DO l=1,llm
	711	DO j=1,jjp1
	712	DO i=1,iip1
	713	q(i,j,l,iq)=q(i,j,l,iq)*patm/patmn
	714	ENDDO
	715	ENDDO
	716	ENDDO
	717	ENDDO
	718	else
	719	write(,) 'OK : conservation of tracer mixing ratio'
	720	end if
	721
	722	c Correction pour la pression au niveau de la mer
	723	yes=' '
	724	do while ((yes.ne.'y').and.(yes.ne.'n'))
	725	write(,) 'Do you wish fix pressure at sea level (y)',
	726	& ' or not (n) ?'
	727	read(*,fmt='(a)') yes
	728	end do
	729
	730	if (yes.eq.'y') then
	731	write(,) 'Value?'
	732	read(,,iostat=ierr) psea
	733	DO i=1,iip1
	734	DO j=1,jjp1
	735	ps(i,j)=psea
	736
	737	ENDDO
	738	ENDDO
	739	write(,) 'psea=',psea
	740	else
	741	write(,) 'no'
	742	end if
	743
	744
	745	c emis : change surface emissivity (added by RW)
	746	c ----------------------------------------------
	747	else if (trim(modif).eq.'emis') then
	748
	749	print*,'new value?'
	750	read(,) emisread
	751
	752	do i=1,ngridmx
	753	emis(i)=emisread
	754	enddo
	755
	756	c qname : change tracer name
	757	c --------------------------
	758	else if (trim(modif).eq.'qname') then
	759	yes='y'
	760	do while (yes.eq.'y')
	761	write(,) 'Which tracer name do you want to change ?'
	762	do iq=1,nqmx
	763	write(*,'(i3,a3,a20)')iq,' : ',trim(tnom(iq))
	764	enddo
	765	write(*,'(a35,i3)')
	766	& '(enter tracer number; between 1 and ',nqmx
	767	write(,)' or any other value to quit this option)'
	768	read(,) iq
	769	if ((iq.ge.1).and.(iq.le.nqmx)) then
	770	write(,)'Change tracer name ',trim(tnom(iq)),' to ?'
	771	read(,) txt
	772	tnom(iq)=txt
	773	write(,)'Do you want to change another tracer name (y/n)?'
	774	read(*,'(a)') yes
	775	else
	776	! inapropiate value of iq; quit this option
	777	yes='n'
	778	endif ! of if ((iq.ge.1).and.(iq.le.nqmx))
	779	enddo ! of do while (yes.ne.'y')
	780
	781	c q=0 : set tracers to zero
	782	c -------------------------
[699]	783	else if (trim(modif).eq.'q=0') then
[135]	784	c mise a 0 des q (traceurs)
	785	write(,) 'Tracers set to 0 (1.E-30 in fact)'
	786	DO iq =1, nqmx
	787	DO l=1,llm
	788	DO j=1,jjp1
	789	DO i=1,iip1
	790	q(i,j,l,iq)=1.e-30
	791	ENDDO
	792	ENDDO
	793	ENDDO
	794	ENDDO
	795
	796	c set surface tracers to zero
	797	DO iq =1, nqmx
	798	DO ig=1,ngridmx
	799	qsurf(ig,iq)=0.
	800	ENDDO
	801	ENDDO
	802
	803	c q=x : initialise tracer manually
	804	c --------------------------------
[699]	805	else if (trim(modif).eq.'q=x') then
[135]	806	write(,) 'Which tracer do you want to modify ?'
	807	do iq=1,nqmx
	808	write(,)iq,' : ',trim(tnom(iq))
	809	enddo
	810	write(,) '(choose between 1 and ',nqmx,')'
	811	read(,) iq
	812	write(,)'mixing ratio of tracer ',trim(tnom(iq)),
	813	& ' ? (kg/kg)'
	814	read(,) val
	815	DO l=1,llm
	816	DO j=1,jjp1
	817	DO i=1,iip1
	818	q(i,j,l,iq)=val
	819	ENDDO
	820	ENDDO
	821	ENDDO
	822	write(,) 'SURFACE value of tracer ',trim(tnom(iq)),
	823	& ' ? (kg/m2)'
	824	read(,) val
	825	DO ig=1,ngridmx
	826	qsurf(ig,iq)=val
	827	ENDDO
	828
[699]	829	c q=profile : initialize tracer with a given profile
	830	c --------------------------------------------------
	831	else if (trim(modif) .eq. 'q=profile') then
	832	write(,) 'Tracer profile will be sought in ASCII file'
	833	write(,) "'profile_tracer' where 'tracer' is tracer name"
	834	write(,) "(one value per line in file; starting with"
	835	write(,) "surface value, the 1st atmospheric layer"
	836	write(,) "followed by 2nd, etc. up to top of atmosphere)"
	837	write(,) 'Which tracer do you want to set?'
	838	do iq=1,nqmx
	839	write(,)iq,' : ',trim(tnom(iq))
	840	enddo
	841	write(,) '(choose between 1 and ',nqmx,')'
	842	read(,) iq
	843	if ((iq.lt.1).or.(iq.gt.nqmx)) then
	844	! wrong value for iq, go back to menu
	845	write(,) "wrong input value:",iq
	846	cycle
	847	endif
	848	! look for input file 'profile_tracer'
	849	txt="profile_"//trim(tnom(iq))
	850	open(41,file=trim(txt),status='old',form='formatted',
	851	& iostat=ierr)
	852	if (ierr.eq.0) then
	853	! OK, found file 'profile_...', load the profile
	854	do l=1,llm+1
	855	read(41,*,iostat=ierr) profile(l)
	856	if (ierr.ne.0) then ! something went wrong
	857	exit ! quit loop
	858	endif
	859	enddo
	860	if (ierr.eq.0) then
	861	! initialize tracer values
	862	qsurf(:,iq)=profile(1)
	863	do l=1,llm
	864	q(:,:,l,iq)=profile(l+1)
	865	enddo
	866	write(,)'OK, tracer ',trim(tnom(iq)),' initialized ',
	867	& 'using values from file ',trim(txt)
	868	else
	869	write(,)'problem reading file ',trim(txt),' !'
	870	write(,)'No modifications to tracer ',trim(tnom(iq))
	871	endif
	872	else
	873	write(,)'Could not find file ',trim(txt),' !'
	874	write(,)'No modifications to tracer ',trim(tnom(iq))
	875	endif
	876
[135]	877
	878	c wetstart : wet atmosphere with a north to south gradient
	879	c --------------------------------------------------------
[699]	880	else if (trim(modif) .eq. 'wetstart') then
[135]	881	! check that there is indeed a water vapor tracer
	882	if (igcm_h2o_vap.eq.0) then
	883	write(,) "No water vapour tracer! Can't use this option"
	884	stop
	885	endif
	886	DO l=1,llm
	887	DO j=1,jjp1
	888	DO i=1,iip1
	889	q(i,j,l,igcm_h2o_vap)=150.e-6 * (rlatu(j)+pi/2.) / pi
	890	ENDDO
	891	ENDDO
	892	ENDDO
	893
	894	write(,) 'Water mass mixing ratio at north pole='
	895	* ,q(1,1,1,igcm_h2o_vap)
	896	write(,) '---------------------------south pole='
	897	* ,q(1,jjp1,1,igcm_h2o_vap)
	898
	899	c noglacier : remove tropical water ice (to initialize high res sim)
	900	c --------------------------------------------------
[699]	901	else if (trim(modif) .eq. 'noglacier') then
[253]	902	if (igcm_h2o_ice.eq.0) then
	903	write(,) "No water ice tracer! Can't use this option"
	904	stop
	905	endif
[135]	906	do ig=1,ngridmx
	907	j=(ig-2)/iim +2
	908	if(ig.eq.1) j=1
	909	write(,) 'OK: remove surface ice for \|lat\|<45'
	910	if (abs(rlatu(j)*180./pi).lt.45.) then
	911	qsurf(ig,igcm_h2o_ice)=0.
	912	end if
	913	end do
	914
	915
	916	c watercapn : H20 ice on permanent northern cap
	917	c --------------------------------------------------
[699]	918	else if (trim(modif) .eq. 'watercapn') then
[253]	919	if (igcm_h2o_ice.eq.0) then
	920	write(,) "No water ice tracer! Can't use this option"
	921	stop
	922	endif
[135]	923
[253]	924	DO j=1,jjp1
	925	DO i=1,iim
	926	ig=1+(j-2)*iim +i
	927	if(j.eq.1) ig=1
	928	if(j.eq.jjp1) ig=ngridmx
	929
	930	if (rlatu(j)*180./pi.gt.80.) then
[367]	931	qsurf(ig,igcm_h2o_ice)=3.4e3
	932	!do isoil=1,nsoilmx
	933	! ith(i,j,isoil) = 18000. ! thermal inertia
	934	!enddo
[253]	935	write(,)' ==> Ice mesh North boundary (deg)= ',
	936	& rlatv(j-1)*180./pi
	937	end if
	938	ENDDO
	939	ENDDO
	940	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
	941
	942	c$$$ do ig=1,ngridmx
	943	c$$$ j=(ig-2)/iim +2
	944	c$$$ if(ig.eq.1) j=1
	945	c$$$ if (rlatu(j)*180./pi.gt.80.) then
	946	c$$$
	947	c$$$ qsurf(ig,igcm_h2o_ice)=1.e5
	948	c$$$ qsurf(ig,igcm_h2o_vap)=0.0!1.e5
	949	c$$$
	950	c$$$ write(,) 'ig=',ig,' H2O ice mass (kg/m2)= ',
	951	c$$$ & qsurf(ig,igcm_h2o_ice)
	952	c$$$
	953	c$$$ write(,)' ==> Ice mesh South boundary (deg)= ',
	954	c$$$ & rlatv(j)*180./pi
	955	c$$$ end if
	956	c$$$ enddo
	957
[135]	958	c watercaps : H20 ice on permanent southern cap
	959	c -------------------------------------------------
[699]	960	else if (trim(modif) .eq. 'watercaps') then
[253]	961	if (igcm_h2o_ice.eq.0) then
	962	write(,) "No water ice tracer! Can't use this option"
	963	stop
	964	endif
[135]	965
[253]	966	DO j=1,jjp1
	967	DO i=1,iim
	968	ig=1+(j-2)*iim +i
	969	if(j.eq.1) ig=1
	970	if(j.eq.jjp1) ig=ngridmx
	971
	972	if (rlatu(j)*180./pi.lt.-80.) then
[367]	973	qsurf(ig,igcm_h2o_ice)=3.4e3
	974	!do isoil=1,nsoilmx
	975	! ith(i,j,isoil) = 18000. ! thermal inertia
	976	!enddo
[253]	977	write(,)' ==> Ice mesh South boundary (deg)= ',
	978	& rlatv(j-1)*180./pi
	979	end if
	980	ENDDO
	981	ENDDO
	982	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
	983
	984	c$$$ do ig=1,ngridmx
	985	c$$$ j=(ig-2)/iim +2
	986	c$$$ if(ig.eq.1) j=1
	987	c$$$ if (rlatu(j)*180./pi.lt.-80.) then
	988	c$$$ qsurf(ig,igcm_h2o_ice)=1.e5
	989	c$$$ qsurf(ig,igcm_h2o_vap)=0.0 !1.e5
	990	c$$$
	991	c$$$ write(,) 'ig=',ig,' H2O ice mass (kg/m2)= ',
	992	c$$$ & qsurf(ig,igcm_h2o_ice)
	993	c$$$ write(,)' ==> Ice mesh North boundary (deg)= ',
	994	c$$$ & rlatv(j-1)*180./pi
	995	c$$$ end if
	996	c$$$ enddo
	997
	998
	999	c noacglac : H2O ice across highest terrain
[135]	1000	c --------------------------------------------
[699]	1001	else if (trim(modif) .eq. 'noacglac') then
[253]	1002	if (igcm_h2o_ice.eq.0) then
	1003	write(,) "No water ice tracer! Can't use this option"
	1004	stop
	1005	endif
	1006	DO j=1,jjp1
	1007	DO i=1,iim
	1008	ig=1+(j-2)*iim +i
	1009	if(j.eq.1) ig=1
	1010	if(j.eq.jjp1) ig=ngridmx
	1011
	1012	if(phis(i,j).gt.1000.*g)then
	1013	alb(i,j) = 0.5 ! snow value
	1014	do isoil=1,nsoilmx
	1015	ith(i,j,isoil) = 18000. ! thermal inertia
	1016	! this leads to rnat set to 'ocean' in physiq.F90
	1017	! actually no, because it is soil not surface
	1018	enddo
	1019	endif
	1020	ENDDO
	1021	ENDDO
	1022	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
	1023	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
	1024	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
	1025
	1026
	1027
	1028	c oborealis : H2O oceans across Vastitas Borealis
	1029	c -----------------------------------------------
[699]	1030	else if (trim(modif) .eq. 'oborealis') then
[253]	1031	if (igcm_h2o_ice.eq.0) then
	1032	write(,) "No water ice tracer! Can't use this option"
	1033	stop
	1034	endif
[135]	1035	DO j=1,jjp1
	1036	DO i=1,iim
	1037	ig=1+(j-2)*iim +i
	1038	if(j.eq.1) ig=1
	1039	if(j.eq.jjp1) ig=ngridmx
	1040
[253]	1041	if(phis(i,j).lt.-4000.*g)then
	1042	! if( (phis(i,j).lt.-4000.*g)
	1043	! & .and. (rlatu(j)*180./pi.lt.0.) )then ! south hemisphere only
	1044	! phis(i,j)=-8000.0*g ! proper ocean
	1045	phis(i,j)=-4000.0*g ! scanty ocean
	1046
	1047	alb(i,j) = 0.07 ! oceanic value
	1048	do isoil=1,nsoilmx
	1049	ith(i,j,isoil) = 18000. ! thermal inertia
	1050	! this leads to rnat set to 'ocean' in physiq.F90
	1051	! actually no, because it is soil not surface
	1052	enddo
[135]	1053	endif
	1054	ENDDO
	1055	ENDDO
[253]	1056	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
	1057	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
	1058	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
[135]	1059
[253]	1060	c iborealis : H2O ice in Northern plains
	1061	c --------------------------------------
[699]	1062	else if (trim(modif) .eq. 'iborealis') then
[253]	1063	if (igcm_h2o_ice.eq.0) then
	1064	write(,) "No water ice tracer! Can't use this option"
	1065	stop
	1066	endif
	1067	DO j=1,jjp1
	1068	DO i=1,iim
	1069	ig=1+(j-2)*iim +i
	1070	if(j.eq.1) ig=1
	1071	if(j.eq.jjp1) ig=ngridmx
	1072
	1073	if(phis(i,j).lt.-4000.*g)then
[367]	1074	!qsurf(ig,igcm_h2o_ice)=1.e3
	1075	qsurf(ig,igcm_h2o_ice)=241.4 ! to make total 33 kg m^-2
[253]	1076	endif
	1077	ENDDO
	1078	ENDDO
	1079	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
	1080	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
	1081	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
	1082
	1083
	1084	c oceanball : H2O liquid everywhere
	1085	c ----------------------------
[699]	1086	else if (trim(modif) .eq. 'oceanball') then
[253]	1087	if (igcm_h2o_ice.eq.0) then
	1088	write(,) "No water ice tracer! Can't use this option"
	1089	stop
	1090	endif
	1091	DO j=1,jjp1
	1092	DO i=1,iim
	1093	ig=1+(j-2)*iim +i
	1094	if(j.eq.1) ig=1
	1095	if(j.eq.jjp1) ig=ngridmx
	1096
	1097	qsurf(ig,igcm_h2o_vap)=0.0 ! for ocean, this is infinite source
	1098	qsurf(ig,igcm_h2o_ice)=0.0
	1099	alb(i,j) = 0.07 ! ocean value
	1100
	1101	do isoil=1,nsoilmx
	1102	ith(i,j,isoil) = 18000. ! thermal inertia
	1103	!ith(i,j,isoil) = 50. ! extremely low for test
	1104	! this leads to rnat set to 'ocean' in physiq.F90
	1105	enddo
	1106
	1107	ENDDO
	1108	ENDDO
	1109	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
	1110	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
	1111	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
	1112
	1113	c iceball : H2O ice everywhere
	1114	c ----------------------------
[699]	1115	else if (trim(modif) .eq. 'iceball') then
[253]	1116	if (igcm_h2o_ice.eq.0) then
	1117	write(,) "No water ice tracer! Can't use this option"
	1118	stop
	1119	endif
	1120	DO j=1,jjp1
	1121	DO i=1,iim
	1122	ig=1+(j-2)*iim +i
	1123	if(j.eq.1) ig=1
	1124	if(j.eq.jjp1) ig=ngridmx
	1125
	1126	qsurf(ig,igcm_h2o_vap)=-50. ! for ocean, this is infinite source
	1127	qsurf(ig,igcm_h2o_ice)=50. ! == to 0.05 m of oceanic ice
	1128	alb(i,j) = 0.6 ! ice albedo value
	1129
	1130	do isoil=1,nsoilmx
	1131	!ith(i,j,isoil) = 18000. ! thermal inertia
	1132	! this leads to rnat set to 'ocean' in physiq.F90
	1133	enddo
	1134
	1135	ENDDO
	1136	ENDDO
	1137	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
	1138	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
	1139
[135]	1140	c isotherm : Isothermal temperatures and no winds
[253]	1141	c -----------------------------------------------
[699]	1142	else if (trim(modif) .eq. 'isotherm') then
[135]	1143
	1144	write(,)'Isothermal temperature of the atmosphere,
	1145	& surface and subsurface'
	1146	write(,) 'Value of this temperature ? :'
	1147	203 read(,,iostat=ierr) Tiso
	1148	if(ierr.ne.0) goto 203
	1149
[699]	1150	tsurf(1:ngridmx)=Tiso
	1151
	1152	tsoil(1:ngridmx,1:nsoilmx)=Tiso
	1153
	1154	Tset(1:iip1,1:jjp1,1:llm)=Tiso
[135]	1155	flagtset=.true.
[699]	1156
	1157	ucov(1:iip1,1:jjp1,1:llm)=0
	1158	vcov(1:iip1,1:jjm,1:llm)=0
	1159	q2(1:ngridmx,1:nlayermx+1)=0
[135]	1160
[253]	1161	c radequi : Radiative equilibrium profile of temperatures and no winds
	1162	c --------------------------------------------------------------------
[699]	1163	else if (trim(modif) .eq. 'radequi') then
[135]	1164
[253]	1165	write(,)'radiative equilibrium temperature profile'
	1166
	1167	do ig=1, ngridmx
	1168	teque= 335.0-60.0sin(latfi(ig))sin(latfi(ig))-
	1169	& 10.0cos(latfi(ig))cos(latfi(ig))
	1170	tsurf(ig) = MAX(220.0,teque)
	1171	end do
	1172	do l=2,nsoilmx
	1173	do ig=1, ngridmx
	1174	tsoil(ig,l) = tsurf(ig)
	1175	end do
	1176	end do
	1177	DO j=1,jjp1
	1178	DO i=1,iim
	1179	Do l=1,llm
	1180	teque=335.-60.0sin(rlatu(j))sin(rlatu(j))-
	1181	& 10.0cos(rlatu(j))cos(rlatu(j))
	1182	Tset(i,j,l)=MAX(220.0,teque)
	1183	end do
	1184	end do
	1185	end do
	1186	flagtset=.true.
[699]	1187	ucov(1:iip1,1:jjp1,1:llm)=0
	1188	vcov(1:iip1,1:jjm,1:llm)=0
	1189	q2(1:ngridmx,1:nlayermx+1)=0
[253]	1190
[135]	1191	c coldstart : T set 1K above CO2 frost point and no winds
	1192	c ------------------------------------------------
[699]	1193	else if (trim(modif) .eq. 'coldstart') then
[135]	1194
	1195	write(,)'set temperature of the atmosphere,'
	1196	&,'surface and subsurface how many degrees above CO2 frost point?'
	1197	204 read(,,iostat=ierr) Tabove
	1198	if(ierr.ne.0) goto 204
	1199
	1200	DO j=1,jjp1
	1201	DO i=1,iim
	1202	ig=1+(j-2)*iim +i
	1203	if(j.eq.1) ig=1
	1204	if(j.eq.jjp1) ig=ngridmx
	1205	tsurf(ig) = (-3167.8)/(log(.01*ps(i,j))-23.23)+Tabove
	1206	END DO
	1207	END DO
	1208	do l=1,nsoilmx
	1209	do ig=1, ngridmx
	1210	tsoil(ig,l) = tsurf(ig)
	1211	end do
	1212	end do
	1213	DO j=1,jjp1
	1214	DO i=1,iim
	1215	Do l=1,llm
	1216	pp = aps(l) +bps(l)*ps(i,j)
	1217	Tset(i,j,l)=(-3167.8)/(log(.01*pp)-23.23)+Tabove
	1218	end do
	1219	end do
	1220	end do
	1221
	1222	flagtset=.true.
[699]	1223	ucov(1:iip1,1:jjp1,1:llm)=0
	1224	vcov(1:iip1,1:jjm,1:llm)=0
	1225	q2(1:ngridmx,1:nlayermx+1)=0
[135]	1226
	1227
	1228	c co2ice=0 : remove CO2 polar ice caps'
	1229	c ------------------------------------------------
[699]	1230	else if (trim(modif) .eq. 'co2ice=0') then
[135]	1231	! check that there is indeed a co2_ice tracer ...
	1232	if (igcm_co2_ice.ne.0) then
	1233	do ig=1,ngridmx
	1234	!co2ice(ig)=0
	1235	qsurf(ig,igcm_co2_ice)=0
	1236	emis(ig)=emis(ngridmx/2)
	1237	end do
	1238	else
	1239	write(,) "Can't remove CO2 ice!! (no co2_ice tracer)"
	1240	endif
	1241
	1242	! therm_ini_s: (re)-set soil thermal inertia to reference surface values
	1243	! ----------------------------------------------------------------------
	1244
[699]	1245	else if (trim(modif).eq.'therm_ini_s') then
[135]	1246	! write(,)"surfithfi(1):",surfithfi(1)
	1247	do isoil=1,nsoilmx
	1248	inertiedat(1:ngridmx,isoil)=surfithfi(1:ngridmx)
	1249	enddo
	1250	write(,)'OK: Soil thermal inertia has been reset to referenc
	1251	&e surface values'
	1252	! write(,)"inertiedat(1,1):",inertiedat(1,1)
	1253	ithfi(:,:)=inertiedat(:,:)
	1254	! recast ithfi() onto ith()
	1255	call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
	1256	! Check:
	1257	! do i=1,iip1
	1258	! do j=1,jjp1
	1259	! do isoil=1,nsoilmx
	1260	! write(77,*) i,j,isoil," ",ith(i,j,isoil)
	1261	! enddo
	1262	! enddo
	1263	! enddo
	1264
	1265
[699]	1266	else
	1267	write(,) ' Unknown (misspelled?) option!!!'
	1268	end if ! of if (trim(modif) .eq. '...') elseif ...
[253]	1269
[135]	1270
	1271
	1272
	1273
	1274
	1275
	1276
	1277
	1278
[253]	1279
[135]	1280	enddo ! of do ! infinite loop on liste of changes
	1281
	1282	999 continue
	1283
	1284
	1285	c=======================================================================
[253]	1286	c Initialisation for cloud fraction and oceanic ice (added by BC 2010)
	1287	c=======================================================================
	1288	DO ig=1, ngridmx
	1289	totalfrac(ig)=0.5
	1290	DO l=1,llm
	1291	cloudfrac(ig,l)=0.5
	1292	ENDDO
[588]	1293	! Ehouarn, march 2012: also add some initialisation for hice
	1294	hice(ig)=0.0
[253]	1295	ENDDO
	1296
	1297	c========================================================
	1298
	1299	! DO ig=1,ngridmx
	1300	! IF(tsurf(ig) .LT. 273.16-1.8) THEN
	1301	! hice(ig)=(273.16-1.8-tsurf(ig))/(273.16-1.8-240)*1.
	1302	! hice(ig)=min(hice(ig),1.0)
	1303	! ENDIF
	1304	! ENDDO
	1305
	1306	c=======================================================================
[135]	1307	c Correct pressure on the new grid (menu 0)
	1308	c=======================================================================
	1309
[588]	1310
[135]	1311	if ((choix_1.eq.0).and.(.not.flagps0)) then
	1312	r = 1000.*8.31/mugaz
	1313
	1314	do j=1,jjp1
	1315	do i=1,iip1
[588]	1316	ps(i,j) = ps(i,j) *
	1317	. exp((phisold_newgrid(i,j)-phis(i,j)) /
	1318	. (t(i,j,1) * r))
[253]	1319	end do
	1320	end do
	1321
[588]	1322	c periodicite de ps en longitude
[253]	1323	do j=1,jjp1
[588]	1324	ps(1,j) = ps(iip1,j)
[135]	1325	end do
	1326	end if
	1327
[588]	1328
[135]	1329	c=======================================================================
	1330	c=======================================================================
	1331
	1332	c=======================================================================
	1333	c Initialisation de la physique / ecriture de newstartfi :
	1334	c=======================================================================
	1335
[588]	1336
[135]	1337	CALL inifilr
	1338	CALL pression(ip1jmp1, ap, bp, ps, p3d)
	1339
	1340	c-----------------------------------------------------------------------
	1341	c Initialisation pks:
	1342	c-----------------------------------------------------------------------
	1343
	1344	CALL exner_hyb(ip1jmp1, ps, p3d, beta, pks, pk, pkf)
	1345	! Calcul de la temperature potentielle teta
	1346
	1347	if (flagtset) then
	1348	DO l=1,llm
	1349	DO j=1,jjp1
	1350	DO i=1,iim
	1351	teta(i,j,l) = Tset(i,j,l) * cpp/pk(i,j,l)
	1352	ENDDO
	1353	teta (iip1,j,l)= teta (1,j,l)
	1354	ENDDO
	1355	ENDDO
	1356	else if (choix_1.eq.0) then
	1357	DO l=1,llm
	1358	DO j=1,jjp1
	1359	DO i=1,iim
	1360	teta(i,j,l) = t(i,j,l) * cpp/pk(i,j,l)
	1361	ENDDO
	1362	teta (iip1,j,l)= teta (1,j,l)
	1363	ENDDO
	1364	ENDDO
	1365	endif
	1366
	1367	C Calcul intermediaire
	1368	c
	1369	if (choix_1.eq.0) then
	1370	CALL massdair( p3d, masse )
	1371	c
	1372	print *,' ALPHAX ',alphax
	1373
	1374	DO l = 1, llm
	1375	DO i = 1, iim
	1376	xppn(i) = aire( i, 1 ) * masse( i , 1 , l )
	1377	xpps(i) = aire( i,jjp1 ) * masse( i , jjp1 , l )
	1378	ENDDO
	1379	xpn = SUM(xppn)/apoln
	1380	xps = SUM(xpps)/apols
	1381	DO i = 1, iip1
	1382	masse( i , 1 , l ) = xpn
	1383	masse( i , jjp1 , l ) = xps
	1384	ENDDO
	1385	ENDDO
	1386	endif
	1387	phis(iip1,:) = phis(1,:)
	1388
	1389	CALL inidissip ( lstardis, nitergdiv, nitergrot, niterh,
	1390	* tetagdiv, tetagrot , tetatemp )
	1391	itau=0
	1392	if (choix_1.eq.0) then
	1393	day_ini=int(date)
	1394	endif
	1395	c
	1396	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
	1397
	1398	CALL caldyn0( itau,ucov,vcov,teta,ps,masse,pk,phis ,
	1399	* phi,w, pbaru,pbarv,day_ini+time )
	1400
[588]	1401
[135]	1402	CALL dynredem0("restart.nc",day_ini,anneeref,phis,nqmx)
	1403	CALL dynredem1("restart.nc",0.0,vcov,ucov,teta,q,nqmx,masse,ps)
	1404	C
	1405	C Ecriture etat initial physique
	1406	C
	1407
[253]	1408	! do ig=1,ngridmx
	1409	! print*,'surface ice in newstart=',qsurf(ig,igcm_h2o_ice)
	1410	! print*,'surface liq in newstart=',qsurf(ig,igcm_h2o_vap)
	1411	! enddo
	1412	! stop
	1413
	1414
[135]	1415	call physdem1("restartfi.nc",lonfi,latfi,nsoilmx,nqmx,
	1416	. dtphys,float(day_ini),
	1417	. time,tsurf,tsoil,emis,q2,qsurf,
[253]	1418	. airefi,albfi,ithfi,zmea,zstd,zsig,zgam,zthe,
	1419	. cloudfrac,totalfrac,hice)
[135]	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

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

Download in other formats:

Original Format