Context Navigation

newstart.F @ 1644

Last change on this file since 1644 was 1644, checked in by jvatant, 8 years ago
Copy the new LMDZ.TITAN from LMDZ.GENERIC + compilation-needed modification (phystd->phytitan ...) + no changes on substance + deftank cleaned
File size: 57.6 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
[1543]	17	use mod_phys_lmdz_para, only: is_parallel, is_sequential,
	18	& is_mpi_root, is_omp_root,
	19	& is_master
[1415]	20	use infotrac, only: infotrac_init, nqtot, tname
[1216]	21	USE tracer_h, ONLY: igcm_co2_ice, igcm_h2o_vap, igcm_h2o_ice
	22	USE comsoil_h, ONLY: nsoilmx, layer, mlayer, inertiedat
	23	USE surfdat_h, ONLY: phisfi, albedodat,
	24	& zmea, zstd, zsig, zgam, zthe
[1470]	25	use datafile_mod, only: datadir, surfdir
[1524]	26	use ioipsl_getin_p_mod, only: getin_p
[1415]	27	use control_mod, only: day_step, iphysiq, anneeref, planet_type
[1216]	28	use phyredem, only: physdem0, physdem1
	29	use iostart, only: open_startphy
[1297]	30	use slab_ice_h, only:noceanmx
[1403]	31	use filtreg_mod, only: inifilr
[1543]	32	USE mod_const_mpi, ONLY: COMM_LMDZ
[1422]	33	USE comvert_mod, ONLY: ap,bp,aps,bps,pa,preff
	34	USE comconst_mod, ONLY: lllm,daysec,dtvr,dtphys,cpp,kappa,
[1543]	35	. rad,omeg,g,r,pi
[1422]	36	USE serre_mod, ONLY: alphax
	37	USE temps_mod, ONLY: day_ini
	38	USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
[1543]	39	use iniphysiq_mod, only: iniphysiq
[135]	40	implicit none
	41
[1543]	42	include "dimensions.h"
[1308]	43	integer, parameter :: ngridmx = (2+(jjm-1)*iim - 1/jjm)
[1543]	44	include "paramet.h"
	45	include "comgeom2.h"
	46	include "comdissnew.h"
	47	include "netcdf.inc"
	48
[135]	49	c=======================================================================
	50	c Declarations
	51	c=======================================================================
	52
	53	c Variables dimension du fichier "start_archive"
	54	c------------------------------------
[1543]	55	CHARACTER relief*3
[135]	56
	57
	58	c Variables pour les lectures NetCDF des fichiers "start_archive"
	59	c--------------------------------------------------
	60	INTEGER nid_dyn, nid_fi,nid,nvarid
	61	INTEGER length
	62	parameter (length = 100)
	63	INTEGER tab0
	64	INTEGER NB_ETATMAX
	65	parameter (NB_ETATMAX = 100)
	66
	67	REAL date
	68	REAL p_rad,p_omeg,p_g,p_cpp,p_mugaz,p_daysec
	69
	70	c Variable histoire
	71	c------------------
	72	REAL vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm) ! vents covariants
	73	REAL phis(iip1,jjp1)
[1216]	74	REAL,ALLOCATABLE :: q(:,:,:,:) ! champs advectes
[135]	75
	76	c autre variables dynamique nouvelle grille
	77	c------------------------------------------
	78	REAL pks(iip1,jjp1)
	79	REAL w(iip1,jjp1,llm+1)
	80	REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
	81	! REAL dv(ip1jm,llm),du(ip1jmp1,llm)
	82	! REAL dh(ip1jmp1,llm),dp(ip1jmp1)
	83	REAL phi(iip1,jjp1,llm)
	84
	85	integer klatdat,klongdat
	86	PARAMETER (klatdat=180,klongdat=360)
	87
	88	c Physique sur grille scalaire
	89	c----------------------------
	90	real zmeaS(iip1,jjp1),zstdS(iip1,jjp1)
	91	real zsigS(iip1,jjp1),zgamS(iip1,jjp1),ztheS(iip1,jjp1)
	92
	93	c variable physique
	94	c------------------
[1543]	95	REAL tsurf(ngridmx) ! surface temperature
	96	REAL,ALLOCATABLE :: tsoil(:,:) ! soil temperature
	97	! REAL co2ice(ngridmx) ! CO2 ice layer
	98	REAL emis(ngridmx) ! surface emissivity
[135]	99	real emisread ! added by RW
[1216]	100	REAL,ALLOCATABLE :: qsurf(:,:)
[1308]	101	REAL q2(ngridmx,llm+1)
[135]	102	! REAL rnaturfi(ngridmx)
	103	real alb(iip1,jjp1),albfi(ngridmx) ! albedos
[1543]	104	real,ALLOCATABLE :: ith(:,:,:),ithfi(:,:) ! thermal inertia (3D)
[135]	105	real surfith(iip1,jjp1),surfithfi(ngridmx) ! surface thermal inertia (2D)
	106	REAL latfi(ngridmx),lonfi(ngridmx),airefi(ngridmx)
	107
	108	INTEGER i,j,l,isoil,ig,idum
	109	real mugaz ! molar mass of the atmosphere
	110
	111	integer ierr
	112
[1297]	113	REAL rnat(ngridmx)
[1543]	114	REAL,ALLOCATABLE :: tslab(:,:) ! slab ocean temperature
[1297]	115	REAL pctsrf_sic(ngridmx) ! sea ice cover
	116	REAL tsea_ice(ngridmx) ! temperature sea_ice
	117	REAL sea_ice(ngridmx) ! mass of sea ice (kg/m2)
	118
[135]	119	c Variables on the new grid along scalar points
	120	c------------------------------------------------------
	121	! REAL p(iip1,jjp1)
	122	REAL t(iip1,jjp1,llm)
	123	REAL tset(iip1,jjp1,llm)
	124	real phisold_newgrid(iip1,jjp1)
	125	REAL :: teta(iip1, jjp1, llm)
	126	REAL :: pk(iip1,jjp1,llm)
	127	REAL :: pkf(iip1,jjp1,llm)
	128	REAL :: ps(iip1, jjp1)
	129	REAL :: masse(iip1,jjp1,llm)
	130	REAL :: xpn,xps,xppn(iim),xpps(iim)
	131	REAL :: p3d(iip1, jjp1, llm+1)
	132	REAL :: beta(iip1,jjp1,llm)
	133	! REAL dteta(ip1jmp1,llm)
	134
	135	c Variable de l'ancienne grille
	136	c------------------------------
	137	real time
	138	real tab_cntrl(100)
	139	real tab_cntrl_bis(100)
	140
	141	c variables diverses
	142	c-------------------
	143	real choix_1,pp
	144	character*80 fichnom
[988]	145	character*250 filestring
[1543]	146	integer Lmodif,iq
[135]	147	character modif*20
	148	real z_reel(iip1,jjp1)
	149	real tsud,albsud,alb_bb,ith_bb,Tiso,Tabove
	150	real ptoto,pcap,patm,airetot,ptotn,patmn,psea
	151	! real ssum
	152	character*1 yes
	153	logical :: flagtset=.false. , flagps0=.false.
[1369]	154	real val, val2, val3, val4 ! to store temporary variables
[135]	155	real :: iceith=2000 ! thermal inertia of subterranean ice
	156
	157	INTEGER :: itau
	158
	159	character(len=20) :: txt ! to store some text
[988]	160	character(len=50) :: surfacefile ! "surface.nc" (or equivalent file)
	161	character(len=150) :: longline
[135]	162	integer :: count
[699]	163	real :: profile(llm+1) ! to store an atmospheric profile + surface value
[135]	164
[253]	165	! added by BC for equilibrium temperature startup
	166	real teque
	167
	168	! added by BC for cloud fraction setup
[1308]	169	REAL hice(ngridmx),cloudfrac(ngridmx,llm)
[253]	170	REAL totalfrac(ngridmx)
	171
[1297]	172
[253]	173	! added by RW for nuketharsis
	174	real fact1
	175	real fact2
	176
	177
[135]	178	c sortie visu pour les champs dynamiques
	179	c---------------------------------------
	180	! INTEGER :: visuid
	181	! real :: time_step,t_ops,t_wrt
	182	! CHARACTER*80 :: visu_file
	183
[535]	184	cpp = 0.
	185	preff = 0.
	186	pa = 0. ! to ensure disaster rather than minor error if we don`t
	187	! make deliberate choice of these values elsewhere.
[135]	188
[1644]	189	planet_type="titan"
[1543]	190
[1216]	191	! initialize "serial/parallel" related stuff
[1543]	192	! (required because we call tabfi() below, before calling iniphysiq)
	193	is_sequential=.true.
	194	is_parallel=.false.
	195	is_mpi_root=.true.
	196	is_omp_root=.true.
	197	is_master=.true.
	198
[1216]	199	! Load tracer number and names:
[1415]	200	call infotrac_init
[1216]	201	! allocate arrays
	202	allocate(q(iip1,jjp1,llm,nqtot))
	203	allocate(qsurf(ngridmx,nqtot))
[1415]	204
[1543]	205	! get value of nsoilmx and allocate corresponding arrays
[1589]	206	! a default value of nsoilmx is set in comsoil_h
[1543]	207	call getin_p("nsoilmx",nsoilmx)
	208
	209	allocate(tsoil(ngridmx,nsoilmx))
	210	allocate(ith(iip1,jjp1,nsoilmx),ithfi(ngridmx,nsoilmx))
	211	allocate(tslab(ngridmx,nsoilmx))
	212
[135]	213	c=======================================================================
	214	c Choice of the start file(s) to use
	215	c=======================================================================
	216	write(,) 'From which kind of files do you want to create new',
	217	. 'start and startfi files'
	218	write(,) ' 0 - from a file start_archive'
	219	write(,) ' 1 - from files start and startfi'
	220
	221	c-----------------------------------------------------------------------
	222	c Open file(s) to modify (start or start_archive)
	223	c-----------------------------------------------------------------------
	224
	225	DO
	226	read(,,iostat=ierr) choix_1
	227	if ((choix_1 /= 0).OR.(choix_1 /=1)) EXIT
	228	ENDDO
	229
	230	c Open start_archive
	231	c ~~~~~~~~~~~~~~~~~~~~~~~~~~
	232	if (choix_1.eq.0) then
	233
	234	write(,) 'Creating start files from:'
	235	write(,) './start_archive.nc'
	236	write(,)
	237	fichnom = 'start_archive.nc'
	238	ierr = NF_OPEN (fichnom, NF_NOWRITE,nid)
	239	IF (ierr.NE.NF_NOERR) THEN
	240	write(6,*)' Problem opening file:',fichnom
	241	write(6,*)' ierr = ', ierr
	242	CALL ABORT
	243	ENDIF
	244	tab0 = 50
	245	Lmodif = 1
	246
	247	c OR open start and startfi files
	248	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	249	else
	250	write(,) 'Creating start files from:'
	251	write(,) './start.nc and ./startfi.nc'
	252	write(,)
	253	fichnom = 'start.nc'
	254	ierr = NF_OPEN (fichnom, NF_NOWRITE,nid_dyn)
	255	IF (ierr.NE.NF_NOERR) THEN
	256	write(6,*)' Problem opening file:',fichnom
	257	write(6,*)' ierr = ', ierr
	258	CALL ABORT
	259	ENDIF
	260
	261	fichnom = 'startfi.nc'
	262	ierr = NF_OPEN (fichnom, NF_NOWRITE,nid_fi)
	263	IF (ierr.NE.NF_NOERR) THEN
	264	write(6,*)' Problem opening file:',fichnom
	265	write(6,*)' ierr = ', ierr
	266	CALL ABORT
	267	ENDIF
	268
	269	tab0 = 0
	270	Lmodif = 0
	271
	272	endif
	273
[649]	274
	275	c=======================================================================
	276	c INITIALISATIONS DIVERSES
	277	c=======================================================================
	278
[787]	279	! Initialize global tracer indexes (stored in tracer.h)
	280	! ... this has to be done before phyetat0
[1216]	281	call initracer(ngridmx,nqtot,tname)
[787]	282
[837]	283	! Take care of arrays in common modules
	284	! ALLOCATE ARRAYS in surfdat_h (if not already done, e.g. when using start_archive)
	285	IF (.not. ALLOCATED(albedodat)) ALLOCATE(albedodat(ngridmx))
	286	IF (.not. ALLOCATED(phisfi)) ALLOCATE(phisfi(ngridmx))
	287	IF (.not. ALLOCATED(zmea)) ALLOCATE(zmea(ngridmx))
	288	IF (.not. ALLOCATED(zstd)) ALLOCATE(zstd(ngridmx))
	289	IF (.not. ALLOCATED(zsig)) ALLOCATE(zsig(ngridmx))
	290	IF (.not. ALLOCATED(zgam)) ALLOCATE(zgam(ngridmx))
	291	IF (.not. ALLOCATED(zthe)) ALLOCATE(zthe(ngridmx))
[787]	292
[135]	293	c-----------------------------------------------------------------------
	294	c Lecture du tableau des parametres du run (pour la dynamique)
	295	c-----------------------------------------------------------------------
	296	if (choix_1.eq.0) then
	297
	298	write(,) 'reading tab_cntrl START_ARCHIVE'
	299	c
	300	ierr = NF_INQ_VARID (nid, "controle", nvarid)
	301	#ifdef NC_DOUBLE
	302	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tab_cntrl)
	303	#else
	304	ierr = NF_GET_VAR_REAL(nid, nvarid, tab_cntrl)
	305	#endif
	306	c
	307	else if (choix_1.eq.1) then
	308
	309	write(,) 'reading tab_cntrl START'
	310	c
	311	ierr = NF_INQ_VARID (nid_dyn, "controle", nvarid)
	312	#ifdef NC_DOUBLE
	313	ierr = NF_GET_VAR_DOUBLE(nid_dyn, nvarid, tab_cntrl)
	314	#else
	315	ierr = NF_GET_VAR_REAL(nid_dyn, nvarid, tab_cntrl)
	316	#endif
	317	c
	318	write(,) 'reading tab_cntrl STARTFI'
	319	c
	320	ierr = NF_INQ_VARID (nid_fi, "controle", nvarid)
	321	#ifdef NC_DOUBLE
	322	ierr = NF_GET_VAR_DOUBLE(nid_fi, nvarid, tab_cntrl_bis)
	323	#else
	324	ierr = NF_GET_VAR_REAL(nid_fi, nvarid, tab_cntrl_bis)
	325	#endif
	326	c
	327	do i=1,50
	328	tab_cntrl(i+50)=tab_cntrl_bis(i)
	329	enddo
[649]	330	write(,) 'printing tab_cntrl', tab_cntrl
	331	do i=1,100
	332	write(,) i,tab_cntrl(i)
	333	enddo
[1543]	334
[649]	335	write(,) 'Reading file START'
	336	fichnom = 'start.nc'
[1421]	337	CALL dynetat0(fichnom,vcov,ucov,teta,q,masse,
[649]	338	. ps,phis,time)
	339
[1589]	340	CALL iniconst
	341	CALL inigeom
	342
	343	! Initialize the physics
	344	CALL iniphysiq(iim,jjm,llm,
	345	& (jjm-1)*iim+2,comm_lmdz,
	346	& daysec,day_ini,dtphys,
	347	& rlatu,rlatv,rlonu,rlonv,
	348	& aire,cu,cv,rad,g,r,cpp,
	349	& 1)
	350
	351	! Lmodif set to 0 to disable modifications possibility in phyeta0
	352	Lmodif=0
[649]	353	write(,) 'Reading file STARTFI'
	354	fichnom = 'startfi.nc'
[1308]	355	CALL phyetat0 (ngridmx,llm,fichnom,tab0,Lmodif,nsoilmx,
	356	. nqtot,day_ini,time,
[649]	357	. tsurf,tsoil,emis,q2,qsurf, !) ! temporary modif by RDW
[1297]	358	. cloudfrac,totalfrac,hice,rnat,pctsrf_sic,tslab,tsea_ice,
	359	. sea_ice)
[649]	360
[905]	361	! copy albedo and soil thermal inertia on (local) physics grid
[649]	362	do i=1,ngridmx
	363	albfi(i) = albedodat(i)
[1543]	364	do j=1,nsoilmx
[649]	365	ithfi(i,j) = inertiedat(i,j)
[1543]	366	enddo
[649]	367	! build a surfithfi(:) using 1st layer of ithfi(:), which might
[905]	368	! be needed later on if reinitializing soil thermal inertia
[649]	369	surfithfi(i)=ithfi(i,1)
	370	enddo
[905]	371	! also copy albedo and soil thermal inertia on (local) dynamics grid
	372	! so that options below can manipulate either (but must then ensure
	373	! to correctly recast things on physics grid)
	374	call gr_fi_dyn(1,ngridmx,iip1,jjp1,albfi,alb)
	375	call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
	376	call gr_fi_dyn(1,ngridmx,iip1,jjp1,surfithfi,surfith)
[135]	377
	378	endif
	379	c-----------------------------------------------------------------------
[1543]	380	c Initialisation des constantes dynamique
[135]	381	c-----------------------------------------------------------------------
	382
	383	kappa = tab_cntrl(9)
	384	etot0 = tab_cntrl(12)
	385	ptot0 = tab_cntrl(13)
	386	ztot0 = tab_cntrl(14)
	387	stot0 = tab_cntrl(15)
	388	ang0 = tab_cntrl(16)
	389	write(,) "Newstart: kappa,etot0,ptot0,ztot0,stot0,ang0"
	390	write(,) kappa,etot0,ptot0,ztot0,stot0,ang0
	391
	392	! for vertical coordinate
	393	preff=tab_cntrl(18) ! reference surface pressure
	394	pa=tab_cntrl(17) ! reference pressure at which coord is purely pressure
	395	!NB: in start_archive files tab_cntrl(17)=tab_cntrl(18)=0
	396	write(,) "Newstart: preff=",preff," pa=",pa
	397	yes=' '
	398	do while ((yes.ne.'y').and.(yes.ne.'n'))
	399	write(,) "Change the values of preff and pa? (y/n)"
	400	read(*,fmt='(a)') yes
	401	end do
	402	if (yes.eq.'y') then
	403	write(,)"New value of reference surface pressure preff?"
	404	write(,)" (for Mars, typically preff=610)"
	405	read(,) preff
	406	write(,)"New value of reference pressure pa for purely"
	407	write(,)"pressure levels (for hybrid coordinates)?"
	408	write(,)" (for Mars, typically pa=20)"
	409	read(,) pa
	410	endif
	411	c-----------------------------------------------------------------------
	412	c Lecture du tab_cntrl et initialisation des constantes physiques
	413	c - pour start: Lmodif = 0 => pas de modifications possibles
	414	c (modif dans le tabfi de readfi + loin)
	415	c - pour start_archive: Lmodif = 1 => modifications possibles
	416	c-----------------------------------------------------------------------
	417	if (choix_1.eq.0) then
[1216]	418	! tabfi requires that input file be first opened by open_startphy(fichnom)
	419	fichnom = 'start_archive.nc'
	420	call open_startphy(fichnom)
[787]	421	call tabfi (ngridmx,nid,Lmodif,tab0,day_ini,lllm,p_rad,
[135]	422	. p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
	423	else if (choix_1.eq.1) then
[1216]	424	fichnom = 'startfi.nc'
	425	call open_startphy(fichnom)
[649]	426	Lmodif=1 ! Lmodif set to 1 to allow modifications in phyeta0
[787]	427	call tabfi (ngridmx,nid_fi,Lmodif,tab0,day_ini,lllm,p_rad,
[135]	428	. p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
	429	endif
	430
[1321]	431	if (p_omeg .eq. -9999.) then
	432	p_omeg = 8.*atan(1.)/p_daysec
	433	print*,"new value of omega",p_omeg
	434	endif
	435
[135]	436	rad = p_rad
	437	omeg = p_omeg
	438	g = p_g
	439	cpp = p_cpp
	440	mugaz = p_mugaz
	441	daysec = p_daysec
	442
[1330]	443	if (p_omeg .eq. -9999.) then
	444	p_omeg = 8.*atan(1.)/p_daysec
	445	print*,"new value of omega",p_omeg
	446	endif
	447
[535]	448	kappa = 8.3141000./(p_mugazp_cpp) ! added by RDW
[135]	449
	450	c=======================================================================
	451	c INITIALISATIONS DIVERSES
	452	c=======================================================================
	453	! Load parameters from run.def file
	454	CALL defrun_new( 99, .TRUE. )
[1589]	455	! Initialize dynamics geometry and co. (which, when using start.nc may
	456	! have changed because of modifications (tabi, preff, pa) above)
[135]	457	CALL iniconst
	458	CALL inigeom
	459	idum=-1
	460	idum=0
	461
[1589]	462	! Initialize the physics for start_archive only
	463	IF (choix_1.eq.0) THEN
[1543]	464	CALL iniphysiq(iim,jjm,llm,
	465	& (jjm-1)*iim+2,comm_lmdz,
	466	& daysec,day_ini,dtphys,
	467	& rlatu,rlatv,rlonu,rlonv,
	468	& aire,cu,cv,rad,g,r,cpp,
	469	& 1)
[1589]	470	ENDIF
[135]	471
	472	c=======================================================================
	473	c lecture topographie, albedo, inertie thermique, relief sous-maille
	474	c=======================================================================
	475
[988]	476	if (choix_1.eq.0) then ! for start_archive files,
	477	! where an external "surface.nc" file is needed
[135]	478
	479	c do while((relief(1:3).ne.'mol').AND.(relief(1:3).ne.'pla'))
	480	c write(,)
	481	c write(,) 'choix du relief (mola,pla)'
	482	c write(,) '(Topographie MGS MOLA, plat)'
	483	c read(*,fmt='(a3)') relief
	484	relief="mola"
	485	c enddo
	486
[988]	487	! First get the correct value of datadir, if not already done:
	488	! default 'datadir' is set in "datafile_mod"
[1315]	489	call getin_p("datadir",datadir)
[988]	490	write(,) 'Available surface data files are:'
[1470]	491	filestring='ls '//trim(datadir)//'/'//
	492	& trim(surfdir)//' \| grep .nc'
[988]	493	call system(filestring)
[1470]	494	! but in ye old days these files were in datadir, so scan it as well
	495	! for the sake of retro-compatibility
	496	filestring='ls '//trim(datadir)//'/'//' \| grep .nc'
	497	call system(filestring)
[135]	498
[988]	499	write(,) ''
	500	write(,) 'Please choose the relevant file',
	501	& ' (e.g. "surface_mars.nc")'
	502	write(,) ' or "none" to not use any (and set planetary'
	503	write(,) ' albedo and surface thermal inertia)'
	504	read(*,fmt='(a50)') surfacefile
[135]	505
[988]	506	if (surfacefile.ne."none") then
[135]	507
[988]	508	CALL datareadnc(relief,surfacefile,phis,alb,surfith,
	509	& zmeaS,zstdS,zsigS,zgamS,ztheS)
	510	else
	511	! specific case when not using a "surface.nc" file
	512	phis(:,:)=0
	513	zmeaS(:,:)=0
	514	zstdS(:,:)=0
	515	zsigS(:,:)=0
	516	zgamS(:,:)=0
	517	ztheS(:,:)=0
	518
	519	write(,) "Enter value of albedo of the bare ground:"
	520	read(,) alb(1,1)
	521	alb(:,:)=alb(1,1)
	522
	523	write(,) "Enter value of thermal inertia of soil:"
	524	read(,) surfith(1,1)
	525	surfith(:,:)=surfith(1,1)
	526
	527	endif ! of if (surfacefile.ne."none")
[135]	528
[988]	529	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
	530	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,surfith,surfithfi)
	531	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
	532	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zmeaS,zmea)
	533	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zstdS,zstd)
	534	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zsigS,zsig)
	535	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zgamS,zgam)
	536	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ztheS,zthe)
	537
	538	endif ! of if (choix_1.eq.0)
	539
	540
[135]	541	c=======================================================================
	542	c Lecture des fichiers (start ou start_archive)
	543	c=======================================================================
	544
	545	if (choix_1.eq.0) then
	546
	547	write(,) 'Reading file START_ARCHIVE'
[1308]	548	CALL lect_start_archive(ngridmx,llm,
	549	& date,tsurf,tsoil,emis,q2,
[1297]	550	& t,ucov,vcov,ps,teta,phisold_newgrid,q,qsurf,
	551	& surfith,nid,
	552	& rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
[135]	553	write(,) "OK, read start_archive file"
[1543]	554	! copy soil thermal inertia
	555	ithfi(:,:)=inertiedat(:,:)
	556
[135]	557	ierr= NF_CLOSE(nid)
	558
[649]	559	else if (choix_1.eq.1) then
	560	!do nothing, start and startfi have already been read
[135]	561	else
	562	CALL exit(1)
	563	endif
	564
	565	dtvr = daysec/FLOAT(day_step)
	566	dtphys = dtvr * FLOAT(iphysiq)
	567
	568	c=======================================================================
	569	c
	570	c=======================================================================
	571
	572	do ! infinite loop on list of changes
	573
	574	write(,)
	575	write(,)
	576	write(,) 'List of possible changes :'
	577	write(,) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
	578	write(,)
	579	write(,) 'flat : no topography ("aquaplanet")'
[1359]	580	write(,) 'set_ps_to_preff : used if changing preff with topo'
[253]	581	write(,) 'nuketharsis : no Tharsis bulge'
[135]	582	write(,) 'bilball : uniform albedo and thermal inertia'
	583	write(,) 'coldspole : cold subsurface and high albedo at S.pole'
	584	write(,) 'qname : change tracer name'
[996]	585	write(,) 't=profile : read temperature profile in profile.in'
[135]	586	write(,) 'q=0 : ALL tracer =zero'
	587	write(,) 'q=x : give a specific uniform value to one tracer'
[699]	588	write(,) 'q=profile : specify a profile for a tracer'
	589	! write(,) 'ini_q : tracers initialisation for chemistry, water an
	590	! $d ice '
	591	! write(,) 'ini_q-H2O : tracers initialisation for chemistry and
	592	! $ice '
	593	! write(,) 'ini_q-iceH2O : tracers initialisation for chemistry on
	594	! $ly '
[135]	595	write(,) 'noglacier : Remove tropical H2O ice if \|lat\|<45'
	596	write(,) 'watercapn : H20 ice on permanent N polar cap '
	597	write(,) 'watercaps : H20 ice on permanent S polar cap '
[253]	598	write(,) 'noacglac : H2O ice across Noachis Terra'
[135]	599	write(,) 'oborealis : H2O ice across Vastitas Borealis'
[253]	600	write(,) 'iceball : Thick ice layer all over surface'
[1369]	601	write(,) 'supercontinent: Create a continent of given Ab and TI'
[135]	602	write(,) 'wetstart : start with a wet atmosphere'
[253]	603	write(,) 'isotherm : Isothermal Temperatures, wind set to zero'
	604	write(,) 'radequi : Earth-like radiative equilibrium temperature
	605	$ profile (lat-alt) and winds set to zero'
[135]	606	write(,) 'coldstart : Start X K above the CO2 frost point and
	607	$set wind to zero (assumes 100% CO2)'
	608	write(,) 'co2ice=0 : remove CO2 polar cap'
	609	write(,) 'ptot : change total pressure'
	610	write(,) 'emis : change surface emissivity'
[253]	611	write(,) 'therm_ini_s : Set soil thermal inertia to reference sur
[135]	612	&face values'
[1297]	613	write(,) 'slab_ocean_0 : initialisation of slab
	614	$ocean variables'
[135]	615
	616	write(,)
	617	write(,) 'Change to perform ?'
	618	write(,) ' (enter keyword or return to end)'
	619	write(,)
	620
	621	read(*,fmt='(a20)') modif
[1216]	622	if (modif(1:1) .eq. ' ')then
	623	exit ! exit loop on changes
	624	endif
[135]	625
	626	write(,)
	627	write(,) trim(modif) , ' : '
	628
	629	c 'flat : no topography ("aquaplanet")'
	630	c -------------------------------------
[699]	631	if (trim(modif) .eq. 'flat') then
[135]	632	c set topo to zero
[699]	633	z_reel(1:iip1,1:jjp1)=0
	634	phis(1:iip1,1:jjp1)=g*z_reel(1:iip1,1:jjp1)
[135]	635	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
	636	write(,) 'topography set to zero.'
	637	write(,) 'WARNING : the subgrid topography parameters',
	638	& ' were not set to zero ! => set calllott to F'
	639
[253]	640	c Choice of surface pressure
[135]	641	yes=' '
	642	do while ((yes.ne.'y').and.(yes.ne.'n'))
	643	write(,) 'Do you wish to choose homogeneous surface',
	644	& 'pressure (y) or let newstart interpolate ',
	645	& ' the previous field (n)?'
	646	read(*,fmt='(a)') yes
	647	end do
	648	if (yes.eq.'y') then
	649	flagps0=.true.
	650	write(,) 'New value for ps (Pa) ?'
	651	201 read(,,iostat=ierr) patm
	652	if(ierr.ne.0) goto 201
[1320]	653	write(,) patm
	654	if (patm.eq.-9999.) then
	655	patm = preff
	656	write(,) "we set patm = preff", preff
	657	endif
[135]	658	write(,)
	659	write(,) ' new ps everywhere (Pa) = ', patm
	660	write(,)
	661	do j=1,jjp1
	662	do i=1,iip1
	663	ps(i,j)=patm
	664	enddo
	665	enddo
	666	end if
	667
[1359]	668	c 'set_ps_to_preff' : used if changing preff with topo
	669	c ----------------------------------------------------
	670	else if (trim(modif) .eq. 'set_ps_to_preff') then
	671	do j=1,jjp1
	672	do i=1,iip1
	673	ps(i,j)=preff
	674	enddo
	675	enddo
	676
[253]	677	c 'nuketharsis : no tharsis bulge for Early Mars'
	678	c ---------------------------------------------
[699]	679	else if (trim(modif) .eq. 'nuketharsis') then
[253]	680
	681	DO j=1,jjp1
	682	DO i=1,iim
	683	ig=1+(j-2)*iim +i
	684	if(j.eq.1) ig=1
	685	if(j.eq.jjp1) ig=ngridmx
	686
	687	fact1=(((rlonv(i)180./pi)+100)*2 +
	688	& (rlatu(j)180./pi)2)/65*2
	689	fact2=exp( -fact1**2.5 )
	690
	691	phis(i,j) = phis(i,j) - (phis(i,j)+4000.g)fact2
	692
	693	! if(phis(i,j).gt.2500.*g)then
	694	! if(rlatu(j)*180./pi.gt.-80.)then ! avoid chopping south polar cap
	695	! phis(i,j)=2500.*g
	696	! endif
	697	! endif
	698
	699	ENDDO
	700	ENDDO
	701	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
	702
	703
	704	c bilball : uniform albedo, thermal inertia
	705	c -----------------------------------------
[699]	706	else if (trim(modif) .eq. 'bilball') then
[135]	707	write(,) 'constante albedo and iner.therm:'
	708	write(,) 'New value for albedo (ex: 0.25) ?'
	709	101 read(,,iostat=ierr) alb_bb
	710	if(ierr.ne.0) goto 101
	711	write(,)
	712	write(,) ' uniform albedo (new value):',alb_bb
	713	write(,)
	714
	715	write(,) 'New value for thermal inertia (eg: 247) ?'
	716	102 read(,,iostat=ierr) ith_bb
	717	if(ierr.ne.0) goto 102
	718	write(,) 'uniform thermal inertia (new value):',ith_bb
	719	DO j=1,jjp1
	720	DO i=1,iip1
[1543]	721	alb(i,j) = alb_bb ! albedo
	722	do isoil=1,nsoilmx
	723	ith(i,j,isoil) = ith_bb ! thermal inertia
	724	enddo
[135]	725	END DO
	726	END DO
	727	! CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ith,ithfi)
	728	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
	729	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
	730
	731	c coldspole : sous-sol de la calotte sud toujours froid
	732	c -----------------------------------------------------
[699]	733	else if (trim(modif) .eq. 'coldspole') then
[135]	734	write(,)'new value for the subsurface temperature',
	735	& ' beneath the permanent southern polar cap ? (eg: 141 K)'
	736	103 read(,,iostat=ierr) tsud
	737	if(ierr.ne.0) goto 103
	738	write(,)
	739	write(,) ' new value of the subsurface temperature:',tsud
	740	c nouvelle temperature sous la calotte permanente
	741	do l=2,nsoilmx
	742	tsoil(ngridmx,l) = tsud
	743	end do
	744
	745
	746	write(,)'new value for the albedo',
	747	& 'of the permanent southern polar cap ? (eg: 0.75)'
	748	104 read(,,iostat=ierr) albsud
	749	if(ierr.ne.0) goto 104
	750	write(,)
	751
	752	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	753	c Option 1: only the albedo of the pole is modified :
	754	albfi(ngridmx)=albsud
	755	write(,) 'ig=',ngridmx,' albedo perennial cap ',
	756	& albfi(ngridmx)
	757
	758	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	759	c Option 2 A haute resolution : coordonnee de la vrai calotte ~
	760	c DO j=1,jjp1
	761	c DO i=1,iip1
	762	c ig=1+(j-2)*iim +i
	763	c if(j.eq.1) ig=1
	764	c if(j.eq.jjp1) ig=ngridmx
	765	c if ((rlatu(j)*180./pi.lt.-84.).and.
	766	c & (rlatu(j)*180./pi.gt.-91.).and.
	767	c & (rlonv(i)*180./pi.gt.-91.).and.
	768	c & (rlonv(i)*180./pi.lt.0.)) then
	769	cc albedo de la calotte permanente fixe a albsud
	770	c alb(i,j)=albsud
	771	c write(,) 'lat=',rlatu(j)*180./pi,
	772	c & ' lon=',rlonv(i)*180./pi
	773	cc fin de la condition sur les limites de la calotte permanente
	774	c end if
	775	c ENDDO
	776	c ENDDO
	777	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	778
	779	c CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
	780
	781
	782	c ptot : Modification of the total pressure: ice + current atmosphere
	783	c -------------------------------------------------------------------
[699]	784	else if (trim(modif).eq.'ptot') then
[135]	785
	786	! check if we have a co2_ice surface tracer:
	787	if (igcm_co2_ice.eq.0) then
	788	write(,) " No surface CO2 ice !"
	789	write(,) " only atmospheric pressure will be considered!"
	790	endif
	791	c calcul de la pression totale glace + atm actuelle
	792	patm=0.
	793	airetot=0.
	794	pcap=0.
	795	DO j=1,jjp1
	796	DO i=1,iim
	797	ig=1+(j-2)*iim +i
	798	if(j.eq.1) ig=1
	799	if(j.eq.jjp1) ig=ngridmx
	800	patm = patm + ps(i,j)*aire(i,j)
	801	airetot= airetot + aire(i,j)
	802	if (igcm_co2_ice.ne.0) then
	803	!pcap = pcap + aire(i,j)co2ice(ig)g
	804	pcap = pcap + aire(i,j)qsurf(ig,igcm_co2_ice)g
	805	endif
	806	ENDDO
	807	ENDDO
	808	ptoto = pcap + patm
	809
	810	print*,'Current total pressure at surface (co2 ice + atm) ',
	811	& ptoto/airetot
	812
	813	print*,'new value?'
	814	read(,) ptotn
	815	ptotn=ptotn*airetot
	816	patmn=ptotn-pcap
	817	print*,'ptoto,patm,ptotn,patmn'
	818	print*,ptoto,patm,ptotn,patmn
	819	print*,'Mult. factor for pressure (atm only)', patmn/patm
	820	do j=1,jjp1
	821	do i=1,iip1
	822	ps(i,j)=ps(i,j)*patmn/patm
	823	enddo
	824	enddo
	825
[253]	826
	827
[135]	828	c Correction pour la conservation des traceurs
	829	yes=' '
	830	do while ((yes.ne.'y').and.(yes.ne.'n'))
	831	write(,) 'Do you wish to conserve tracer total mass (y)',
	832	& ' or tracer mixing ratio (n) ?'
	833	read(*,fmt='(a)') yes
	834	end do
	835
	836	if (yes.eq.'y') then
	837	write(,) 'OK : conservation of tracer total mass'
[1216]	838	DO iq =1, nqtot
[135]	839	DO l=1,llm
	840	DO j=1,jjp1
	841	DO i=1,iip1
	842	q(i,j,l,iq)=q(i,j,l,iq)*patm/patmn
	843	ENDDO
	844	ENDDO
	845	ENDDO
	846	ENDDO
	847	else
	848	write(,) 'OK : conservation of tracer mixing ratio'
	849	end if
	850
	851	c Correction pour la pression au niveau de la mer
	852	yes=' '
	853	do while ((yes.ne.'y').and.(yes.ne.'n'))
	854	write(,) 'Do you wish fix pressure at sea level (y)',
	855	& ' or not (n) ?'
	856	read(*,fmt='(a)') yes
	857	end do
	858
	859	if (yes.eq.'y') then
	860	write(,) 'Value?'
[1543]	861	read(,,iostat=ierr) psea
[135]	862	DO i=1,iip1
	863	DO j=1,jjp1
	864	ps(i,j)=psea
	865
	866	ENDDO
	867	ENDDO
[1543]	868	write(,) 'psea=',psea
[135]	869	else
	870	write(,) 'no'
	871	end if
	872
	873
	874	c emis : change surface emissivity (added by RW)
	875	c ----------------------------------------------
	876	else if (trim(modif).eq.'emis') then
	877
	878	print*,'new value?'
	879	read(,) emisread
	880
	881	do i=1,ngridmx
	882	emis(i)=emisread
	883	enddo
	884
	885	c qname : change tracer name
	886	c --------------------------
	887	else if (trim(modif).eq.'qname') then
	888	yes='y'
	889	do while (yes.eq.'y')
	890	write(,) 'Which tracer name do you want to change ?'
[1216]	891	do iq=1,nqtot
	892	write(*,'(i3,a3,a20)')iq,' : ',trim(tname(iq))
[135]	893	enddo
	894	write(*,'(a35,i3)')
[1216]	895	& '(enter tracer number; between 1 and ',nqtot
[135]	896	write(,)' or any other value to quit this option)'
	897	read(,) iq
[1216]	898	if ((iq.ge.1).and.(iq.le.nqtot)) then
	899	write(,)'Change tracer name ',trim(tname(iq)),' to ?'
[135]	900	read(,) txt
[1216]	901	tname(iq)=txt
[135]	902	write(,)'Do you want to change another tracer name (y/n)?'
	903	read(*,'(a)') yes
	904	else
	905	! inapropiate value of iq; quit this option
	906	yes='n'
[1216]	907	endif ! of if ((iq.ge.1).and.(iq.le.nqtot))
[135]	908	enddo ! of do while (yes.ne.'y')
	909
	910	c q=0 : set tracers to zero
	911	c -------------------------
[699]	912	else if (trim(modif).eq.'q=0') then
[135]	913	c mise a 0 des q (traceurs)
	914	write(,) 'Tracers set to 0 (1.E-30 in fact)'
[1216]	915	DO iq =1, nqtot
[135]	916	DO l=1,llm
	917	DO j=1,jjp1
	918	DO i=1,iip1
	919	q(i,j,l,iq)=1.e-30
	920	ENDDO
	921	ENDDO
	922	ENDDO
	923	ENDDO
	924
	925	c set surface tracers to zero
[1216]	926	DO iq =1, nqtot
[135]	927	DO ig=1,ngridmx
	928	qsurf(ig,iq)=0.
	929	ENDDO
	930	ENDDO
	931
	932	c q=x : initialise tracer manually
	933	c --------------------------------
[699]	934	else if (trim(modif).eq.'q=x') then
[135]	935	write(,) 'Which tracer do you want to modify ?'
[1216]	936	do iq=1,nqtot
	937	write(,)iq,' : ',trim(tname(iq))
[135]	938	enddo
[1216]	939	write(,) '(choose between 1 and ',nqtot,')'
[135]	940	read(,) iq
[1216]	941	write(,)'mixing ratio of tracer ',trim(tname(iq)),
[135]	942	& ' ? (kg/kg)'
	943	read(,) val
	944	DO l=1,llm
	945	DO j=1,jjp1
	946	DO i=1,iip1
	947	q(i,j,l,iq)=val
	948	ENDDO
	949	ENDDO
	950	ENDDO
[1216]	951	write(,) 'SURFACE value of tracer ',trim(tname(iq)),
[135]	952	& ' ? (kg/m2)'
	953	read(,) val
	954	DO ig=1,ngridmx
	955	qsurf(ig,iq)=val
	956	ENDDO
	957
[996]	958	c t=profile : initialize temperature with a given profile
	959	c -------------------------------------------------------
	960	else if (trim(modif) .eq. 't=profile') then
	961	write(,) 'Temperature profile from ASCII file'
	962	write(,) "'profile.in' e.g. 1D output"
	963	write(,) "(one value per line in file; starting with"
	964	write(,) "surface value, the 1st atmospheric layer"
	965	write(,) "followed by 2nd, etc. up to top of atmosphere)"
	966	txt="profile.in"
	967	open(33,file=trim(txt),status='old',form='formatted',
	968	& iostat=ierr)
	969	if (ierr.eq.0) then
	970	! OK, found file 'profile_...', load the profile
	971	do l=1,llm+1
	972	read(33,*,iostat=ierr) profile(l)
	973	write(,) profile(l)
	974	if (ierr.ne.0) then ! something went wrong
	975	exit ! quit loop
	976	endif
	977	enddo
	978	if (ierr.eq.0) then
	979	tsurf(1:ngridmx)=profile(1)
	980	tsoil(1:ngridmx,1:nsoilmx)=profile(1)
	981	do l=1,llm
	982	Tset(1:iip1,1:jjp1,l)=profile(l+1)
	983	flagtset=.true.
	984	enddo
	985	ucov(1:iip1,1:jjp1,1:llm)=0.
	986	vcov(1:iip1,1:jjm,1:llm)=0.
[1308]	987	q2(1:ngridmx,1:llm+1)=0.
[996]	988	else
	989	write(,)'problem reading file ',trim(txt),' !'
	990	write(,)'No modifications to temperature'
	991	endif
	992	else
	993	write(,)'Could not find file ',trim(txt),' !'
	994	write(,)'No modifications to temperature'
	995	endif
	996
[699]	997	c q=profile : initialize tracer with a given profile
	998	c --------------------------------------------------
	999	else if (trim(modif) .eq. 'q=profile') then
	1000	write(,) 'Tracer profile will be sought in ASCII file'
	1001	write(,) "'profile_tracer' where 'tracer' is tracer name"
	1002	write(,) "(one value per line in file; starting with"
	1003	write(,) "surface value, the 1st atmospheric layer"
	1004	write(,) "followed by 2nd, etc. up to top of atmosphere)"
	1005	write(,) 'Which tracer do you want to set?'
[1216]	1006	do iq=1,nqtot
	1007	write(,)iq,' : ',trim(tname(iq))
[699]	1008	enddo
[1216]	1009	write(,) '(choose between 1 and ',nqtot,')'
[699]	1010	read(,) iq
[1216]	1011	if ((iq.lt.1).or.(iq.gt.nqtot)) then
[699]	1012	! wrong value for iq, go back to menu
	1013	write(,) "wrong input value:",iq
	1014	cycle
	1015	endif
	1016	! look for input file 'profile_tracer'
[1216]	1017	txt="profile_"//trim(tname(iq))
[699]	1018	open(41,file=trim(txt),status='old',form='formatted',
	1019	& iostat=ierr)
	1020	if (ierr.eq.0) then
	1021	! OK, found file 'profile_...', load the profile
	1022	do l=1,llm+1
	1023	read(41,*,iostat=ierr) profile(l)
	1024	if (ierr.ne.0) then ! something went wrong
	1025	exit ! quit loop
	1026	endif
	1027	enddo
	1028	if (ierr.eq.0) then
	1029	! initialize tracer values
	1030	qsurf(:,iq)=profile(1)
	1031	do l=1,llm
	1032	q(:,:,l,iq)=profile(l+1)
	1033	enddo
[1216]	1034	write(,)'OK, tracer ',trim(tname(iq)),' initialized '
	1035	& ,'using values from file ',trim(txt)
[699]	1036	else
	1037	write(,)'problem reading file ',trim(txt),' !'
[1216]	1038	write(,)'No modifications to tracer ',trim(tname(iq))
[699]	1039	endif
	1040	else
	1041	write(,)'Could not find file ',trim(txt),' !'
[1216]	1042	write(,)'No modifications to tracer ',trim(tname(iq))
[699]	1043	endif
	1044
[135]	1045
	1046	c wetstart : wet atmosphere with a north to south gradient
	1047	c --------------------------------------------------------
[699]	1048	else if (trim(modif) .eq. 'wetstart') then
[135]	1049	! check that there is indeed a water vapor tracer
	1050	if (igcm_h2o_vap.eq.0) then
	1051	write(,) "No water vapour tracer! Can't use this option"
	1052	stop
	1053	endif
	1054	DO l=1,llm
	1055	DO j=1,jjp1
	1056	DO i=1,iip1
	1057	q(i,j,l,igcm_h2o_vap)=150.e-6 * (rlatu(j)+pi/2.) / pi
	1058	ENDDO
	1059	ENDDO
	1060	ENDDO
	1061
	1062	write(,) 'Water mass mixing ratio at north pole='
	1063	* ,q(1,1,1,igcm_h2o_vap)
	1064	write(,) '---------------------------south pole='
	1065	* ,q(1,jjp1,1,igcm_h2o_vap)
	1066
	1067	c noglacier : remove tropical water ice (to initialize high res sim)
	1068	c --------------------------------------------------
[699]	1069	else if (trim(modif) .eq. 'noglacier') then
[253]	1070	if (igcm_h2o_ice.eq.0) then
	1071	write(,) "No water ice tracer! Can't use this option"
	1072	stop
	1073	endif
[135]	1074	do ig=1,ngridmx
	1075	j=(ig-2)/iim +2
	1076	if(ig.eq.1) j=1
	1077	write(,) 'OK: remove surface ice for \|lat\|<45'
	1078	if (abs(rlatu(j)*180./pi).lt.45.) then
	1079	qsurf(ig,igcm_h2o_ice)=0.
	1080	end if
	1081	end do
	1082
	1083
	1084	c watercapn : H20 ice on permanent northern cap
	1085	c --------------------------------------------------
[699]	1086	else if (trim(modif) .eq. 'watercapn') 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
[135]	1091
[253]	1092	DO j=1,jjp1
	1093	DO i=1,iim
	1094	ig=1+(j-2)*iim +i
	1095	if(j.eq.1) ig=1
	1096	if(j.eq.jjp1) ig=ngridmx
	1097
	1098	if (rlatu(j)*180./pi.gt.80.) then
[367]	1099	qsurf(ig,igcm_h2o_ice)=3.4e3
	1100	!do isoil=1,nsoilmx
	1101	! ith(i,j,isoil) = 18000. ! thermal inertia
	1102	!enddo
[253]	1103	write(,)' ==> Ice mesh North boundary (deg)= ',
	1104	& rlatv(j-1)*180./pi
	1105	end if
	1106	ENDDO
	1107	ENDDO
	1108	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
	1109
	1110	c$$$ do ig=1,ngridmx
	1111	c$$$ j=(ig-2)/iim +2
	1112	c$$$ if(ig.eq.1) j=1
	1113	c$$$ if (rlatu(j)*180./pi.gt.80.) then
	1114	c$$$
	1115	c$$$ qsurf(ig,igcm_h2o_ice)=1.e5
	1116	c$$$ qsurf(ig,igcm_h2o_vap)=0.0!1.e5
	1117	c$$$
	1118	c$$$ write(,) 'ig=',ig,' H2O ice mass (kg/m2)= ',
	1119	c$$$ & qsurf(ig,igcm_h2o_ice)
	1120	c$$$
	1121	c$$$ write(,)' ==> Ice mesh South boundary (deg)= ',
	1122	c$$$ & rlatv(j)*180./pi
	1123	c$$$ end if
	1124	c$$$ enddo
	1125
[135]	1126	c watercaps : H20 ice on permanent southern cap
	1127	c -------------------------------------------------
[699]	1128	else if (trim(modif) .eq. 'watercaps') then
[253]	1129	if (igcm_h2o_ice.eq.0) then
	1130	write(,) "No water ice tracer! Can't use this option"
	1131	stop
	1132	endif
[135]	1133
[253]	1134	DO j=1,jjp1
	1135	DO i=1,iim
	1136	ig=1+(j-2)*iim +i
	1137	if(j.eq.1) ig=1
	1138	if(j.eq.jjp1) ig=ngridmx
	1139
	1140	if (rlatu(j)*180./pi.lt.-80.) then
[367]	1141	qsurf(ig,igcm_h2o_ice)=3.4e3
	1142	!do isoil=1,nsoilmx
	1143	! ith(i,j,isoil) = 18000. ! thermal inertia
	1144	!enddo
[253]	1145	write(,)' ==> Ice mesh South boundary (deg)= ',
	1146	& rlatv(j-1)*180./pi
	1147	end if
	1148	ENDDO
	1149	ENDDO
	1150	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
	1151
	1152	c$$$ do ig=1,ngridmx
	1153	c$$$ j=(ig-2)/iim +2
	1154	c$$$ if(ig.eq.1) j=1
	1155	c$$$ if (rlatu(j)*180./pi.lt.-80.) then
	1156	c$$$ qsurf(ig,igcm_h2o_ice)=1.e5
	1157	c$$$ qsurf(ig,igcm_h2o_vap)=0.0 !1.e5
	1158	c$$$
	1159	c$$$ write(,) 'ig=',ig,' H2O ice mass (kg/m2)= ',
	1160	c$$$ & qsurf(ig,igcm_h2o_ice)
	1161	c$$$ write(,)' ==> Ice mesh North boundary (deg)= ',
	1162	c$$$ & rlatv(j-1)*180./pi
	1163	c$$$ end if
	1164	c$$$ enddo
	1165
	1166
	1167	c noacglac : H2O ice across highest terrain
[135]	1168	c --------------------------------------------
[699]	1169	else if (trim(modif) .eq. 'noacglac') then
[253]	1170	if (igcm_h2o_ice.eq.0) then
	1171	write(,) "No water ice tracer! Can't use this option"
	1172	stop
	1173	endif
	1174	DO j=1,jjp1
	1175	DO i=1,iim
	1176	ig=1+(j-2)*iim +i
	1177	if(j.eq.1) ig=1
	1178	if(j.eq.jjp1) ig=ngridmx
	1179
	1180	if(phis(i,j).gt.1000.*g)then
	1181	alb(i,j) = 0.5 ! snow value
	1182	do isoil=1,nsoilmx
	1183	ith(i,j,isoil) = 18000. ! thermal inertia
	1184	! this leads to rnat set to 'ocean' in physiq.F90
	1185	! actually no, because it is soil not surface
	1186	enddo
	1187	endif
	1188	ENDDO
	1189	ENDDO
	1190	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
	1191	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
	1192	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
	1193
	1194
	1195
	1196	c oborealis : H2O oceans across Vastitas Borealis
	1197	c -----------------------------------------------
[699]	1198	else if (trim(modif) .eq. 'oborealis') then
[253]	1199	if (igcm_h2o_ice.eq.0) then
	1200	write(,) "No water ice tracer! Can't use this option"
	1201	stop
	1202	endif
[135]	1203	DO j=1,jjp1
	1204	DO i=1,iim
	1205	ig=1+(j-2)*iim +i
	1206	if(j.eq.1) ig=1
	1207	if(j.eq.jjp1) ig=ngridmx
	1208
[253]	1209	if(phis(i,j).lt.-4000.*g)then
	1210	! if( (phis(i,j).lt.-4000.*g)
	1211	! & .and. (rlatu(j)*180./pi.lt.0.) )then ! south hemisphere only
	1212	! phis(i,j)=-8000.0*g ! proper ocean
	1213	phis(i,j)=-4000.0*g ! scanty ocean
	1214
	1215	alb(i,j) = 0.07 ! oceanic value
	1216	do isoil=1,nsoilmx
	1217	ith(i,j,isoil) = 18000. ! thermal inertia
	1218	! this leads to rnat set to 'ocean' in physiq.F90
	1219	! actually no, because it is soil not surface
	1220	enddo
[135]	1221	endif
	1222	ENDDO
	1223	ENDDO
[253]	1224	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
	1225	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
	1226	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
[135]	1227
[253]	1228	c iborealis : H2O ice in Northern plains
	1229	c --------------------------------------
[699]	1230	else if (trim(modif) .eq. 'iborealis') then
[253]	1231	if (igcm_h2o_ice.eq.0) then
	1232	write(,) "No water ice tracer! Can't use this option"
	1233	stop
	1234	endif
	1235	DO j=1,jjp1
	1236	DO i=1,iim
	1237	ig=1+(j-2)*iim +i
	1238	if(j.eq.1) ig=1
	1239	if(j.eq.jjp1) ig=ngridmx
	1240
	1241	if(phis(i,j).lt.-4000.*g)then
[367]	1242	!qsurf(ig,igcm_h2o_ice)=1.e3
	1243	qsurf(ig,igcm_h2o_ice)=241.4 ! to make total 33 kg m^-2
[253]	1244	endif
	1245	ENDDO
	1246	ENDDO
	1247	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
	1248	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
	1249	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
	1250
	1251
	1252	c oceanball : H2O liquid everywhere
	1253	c ----------------------------
[699]	1254	else if (trim(modif) .eq. 'oceanball') then
[253]	1255	if (igcm_h2o_ice.eq.0) then
	1256	write(,) "No water ice tracer! Can't use this option"
	1257	stop
	1258	endif
	1259	DO j=1,jjp1
	1260	DO i=1,iim
	1261	ig=1+(j-2)*iim +i
	1262	if(j.eq.1) ig=1
	1263	if(j.eq.jjp1) ig=ngridmx
	1264
	1265	qsurf(ig,igcm_h2o_vap)=0.0 ! for ocean, this is infinite source
	1266	qsurf(ig,igcm_h2o_ice)=0.0
	1267	alb(i,j) = 0.07 ! ocean value
	1268
	1269	do isoil=1,nsoilmx
	1270	ith(i,j,isoil) = 18000. ! thermal inertia
	1271	!ith(i,j,isoil) = 50. ! extremely low for test
	1272	! this leads to rnat set to 'ocean' in physiq.F90
	1273	enddo
	1274
	1275	ENDDO
	1276	ENDDO
	1277	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
	1278	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
	1279	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
	1280
	1281	c iceball : H2O ice everywhere
	1282	c ----------------------------
[699]	1283	else if (trim(modif) .eq. 'iceball') then
[253]	1284	if (igcm_h2o_ice.eq.0) then
	1285	write(,) "No water ice tracer! Can't use this option"
	1286	stop
	1287	endif
	1288	DO j=1,jjp1
	1289	DO i=1,iim
	1290	ig=1+(j-2)*iim +i
	1291	if(j.eq.1) ig=1
	1292	if(j.eq.jjp1) ig=ngridmx
	1293
	1294	qsurf(ig,igcm_h2o_vap)=-50. ! for ocean, this is infinite source
	1295	qsurf(ig,igcm_h2o_ice)=50. ! == to 0.05 m of oceanic ice
	1296	alb(i,j) = 0.6 ! ice albedo value
	1297
	1298	do isoil=1,nsoilmx
	1299	!ith(i,j,isoil) = 18000. ! thermal inertia
	1300	! this leads to rnat set to 'ocean' in physiq.F90
	1301	enddo
	1302
	1303	ENDDO
	1304	ENDDO
	1305	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
	1306	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
	1307
[1369]	1308	c supercontinent : H2O ice everywhere
	1309	c ----------------------------
	1310	else if (trim(modif) .eq. 'supercontinent') then
	1311	write(,) 'Minimum longitude (-180,180)'
	1312	read(,) val
	1313	write(,) 'Maximum longitude (-180,180)'
	1314	read(,) val2
	1315	write(,) 'Minimum latitude (-90,90)'
	1316	read(,) val3
	1317	write(,) 'Maximum latitude (-90,90)'
	1318	read(,) val4
	1319
	1320	do j=1,jjp1
	1321	do i=1,iip1
	1322	ig=1+(j-2)*iim +i
	1323	if(j.eq.1) ig=1
	1324	if(j.eq.jjp1) ig=ngridmx
	1325
	1326	c Supercontinent:
[1370]	1327	if (((rlatu(j)*180./pi.le.val4).and.
	1328	& (rlatu(j)*180./pi.ge.val3).and.
	1329	& (rlonv(i)*180./pi.le.val2).and.
	1330	& (rlonv(i)*180./pi.ge.val))) then
[1369]	1331
	1332	rnat(ig)=1.
	1333	alb(i,j) = 0.3
	1334	do isoil=1,nsoilmx
	1335	ith(i,j,isoil) = 2000.
	1336	enddo
	1337	c Ocean:
	1338	else
	1339	rnat(ig)=0.
	1340	alb(i,j) = 0.07
	1341	do isoil=1,nsoilmx
	1342	ith(i,j,isoil) = 18000.
	1343	enddo
	1344	end if
	1345
	1346	enddo
	1347	enddo
	1348	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
	1349	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
	1350
[135]	1351	c isotherm : Isothermal temperatures and no winds
[253]	1352	c -----------------------------------------------
[699]	1353	else if (trim(modif) .eq. 'isotherm') then
[135]	1354
	1355	write(,)'Isothermal temperature of the atmosphere,
	1356	& surface and subsurface'
	1357	write(,) 'Value of this temperature ? :'
	1358	203 read(,,iostat=ierr) Tiso
	1359	if(ierr.ne.0) goto 203
	1360
[699]	1361	tsurf(1:ngridmx)=Tiso
	1362
	1363	tsoil(1:ngridmx,1:nsoilmx)=Tiso
	1364
	1365	Tset(1:iip1,1:jjp1,1:llm)=Tiso
[135]	1366	flagtset=.true.
[1359]	1367
	1368	t(1:iip1,1:jjp1,1:llm)=Tiso
	1369	!! otherwise hydrost. integrations below
	1370	!! use the wrong temperature
	1371	!! -- NB: Tset might be useless
	1372
[699]	1373	ucov(1:iip1,1:jjp1,1:llm)=0
	1374	vcov(1:iip1,1:jjm,1:llm)=0
[1308]	1375	q2(1:ngridmx,1:llm+1)=0
[135]	1376
[253]	1377	c radequi : Radiative equilibrium profile of temperatures and no winds
	1378	c --------------------------------------------------------------------
[699]	1379	else if (trim(modif) .eq. 'radequi') then
[135]	1380
[253]	1381	write(,)'radiative equilibrium temperature profile'
	1382
	1383	do ig=1, ngridmx
	1384	teque= 335.0-60.0sin(latfi(ig))sin(latfi(ig))-
	1385	& 10.0cos(latfi(ig))cos(latfi(ig))
	1386	tsurf(ig) = MAX(220.0,teque)
	1387	end do
	1388	do l=2,nsoilmx
	1389	do ig=1, ngridmx
	1390	tsoil(ig,l) = tsurf(ig)
	1391	end do
	1392	end do
	1393	DO j=1,jjp1
	1394	DO i=1,iim
	1395	Do l=1,llm
	1396	teque=335.-60.0sin(rlatu(j))sin(rlatu(j))-
	1397	& 10.0cos(rlatu(j))cos(rlatu(j))
	1398	Tset(i,j,l)=MAX(220.0,teque)
	1399	end do
	1400	end do
	1401	end do
	1402	flagtset=.true.
[699]	1403	ucov(1:iip1,1:jjp1,1:llm)=0
	1404	vcov(1:iip1,1:jjm,1:llm)=0
[1308]	1405	q2(1:ngridmx,1:llm+1)=0
[253]	1406
[135]	1407	c coldstart : T set 1K above CO2 frost point and no winds
	1408	c ------------------------------------------------
[699]	1409	else if (trim(modif) .eq. 'coldstart') then
[135]	1410
	1411	write(,)'set temperature of the atmosphere,'
	1412	&,'surface and subsurface how many degrees above CO2 frost point?'
	1413	204 read(,,iostat=ierr) Tabove
	1414	if(ierr.ne.0) goto 204
	1415
	1416	DO j=1,jjp1
	1417	DO i=1,iim
	1418	ig=1+(j-2)*iim +i
	1419	if(j.eq.1) ig=1
	1420	if(j.eq.jjp1) ig=ngridmx
	1421	tsurf(ig) = (-3167.8)/(log(.01*ps(i,j))-23.23)+Tabove
	1422	END DO
	1423	END DO
	1424	do l=1,nsoilmx
	1425	do ig=1, ngridmx
	1426	tsoil(ig,l) = tsurf(ig)
	1427	end do
	1428	end do
	1429	DO j=1,jjp1
	1430	DO i=1,iim
	1431	Do l=1,llm
	1432	pp = aps(l) +bps(l)*ps(i,j)
	1433	Tset(i,j,l)=(-3167.8)/(log(.01*pp)-23.23)+Tabove
	1434	end do
	1435	end do
	1436	end do
	1437
	1438	flagtset=.true.
[699]	1439	ucov(1:iip1,1:jjp1,1:llm)=0
	1440	vcov(1:iip1,1:jjm,1:llm)=0
[1308]	1441	q2(1:ngridmx,1:llm+1)=0
[135]	1442
	1443
	1444	c co2ice=0 : remove CO2 polar ice caps'
	1445	c ------------------------------------------------
[699]	1446	else if (trim(modif) .eq. 'co2ice=0') then
[135]	1447	! check that there is indeed a co2_ice tracer ...
	1448	if (igcm_co2_ice.ne.0) then
	1449	do ig=1,ngridmx
	1450	!co2ice(ig)=0
	1451	qsurf(ig,igcm_co2_ice)=0
	1452	emis(ig)=emis(ngridmx/2)
	1453	end do
	1454	else
	1455	write(,) "Can't remove CO2 ice!! (no co2_ice tracer)"
	1456	endif
	1457
	1458	! therm_ini_s: (re)-set soil thermal inertia to reference surface values
	1459	! ----------------------------------------------------------------------
	1460
[1543]	1461	else if (trim(modif) .eq. 'therm_ini_s') then
[135]	1462	! write(,)"surfithfi(1):",surfithfi(1)
[1543]	1463	do isoil=1,nsoilmx
	1464	inertiedat(1:ngridmx,isoil)=surfithfi(1:ngridmx)
	1465	enddo
[135]	1466	write(,)'OK: Soil thermal inertia has been reset to referenc
	1467	&e surface values'
[1543]	1468	! write(,)"inertiedat(1,1):",inertiedat(1,1)
	1469	ithfi(:,:)=inertiedat(:,:)
	1470	! recast ithfi() onto ith()
	1471	call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
[135]	1472	! Check:
	1473	! do i=1,iip1
	1474	! do j=1,jjp1
	1475	! do isoil=1,nsoilmx
	1476	! write(77,*) i,j,isoil," ",ith(i,j,isoil)
	1477	! enddo
	1478	! enddo
[1543]	1479	! enddo
[135]	1480
	1481
	1482
[1297]	1483	c slab_ocean_initialisation
	1484	c ------------------------------------------------
	1485	else if (trim(modif) .eq. 'slab_ocean_0') then
	1486	write(,)'OK: initialisation of slab ocean'
[135]	1487
[1297]	1488	DO ig=1, ngridmx
	1489	rnat(ig)=1.
	1490	tslab(ig,1)=0.
	1491	tslab(ig,2)=0.
	1492	tsea_ice(ig)=0.
	1493	sea_ice(ig)=0.
	1494	pctsrf_sic(ig)=0.
	1495
	1496	if(ithfi(ig,1).GT.10000.)then
	1497	rnat(ig)=0.
	1498	phisfi(ig)=0.
	1499	tsea_ice(ig)=273.15-1.8
	1500	tslab(ig,1)=tsurf(ig)
	1501	tslab(ig,2)=tsurf(ig)!*2./3.+(273.15-1.8)/3.
	1502	endif
[135]	1503
[1297]	1504	ENDDO
	1505	CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
[135]	1506
	1507
	1508
[1297]	1509	else
	1510	write(,) ' Unknown (misspelled?) option!!!'
	1511	end if ! of if (trim(modif) .eq. '...') elseif ...
[135]	1512
	1513
	1514
	1515	enddo ! of do ! infinite loop on liste of changes
	1516
	1517	999 continue
	1518
	1519
	1520	c=======================================================================
[253]	1521	c Initialisation for cloud fraction and oceanic ice (added by BC 2010)
	1522	c=======================================================================
	1523	DO ig=1, ngridmx
	1524	totalfrac(ig)=0.5
	1525	DO l=1,llm
	1526	cloudfrac(ig,l)=0.5
	1527	ENDDO
[588]	1528	! Ehouarn, march 2012: also add some initialisation for hice
	1529	hice(ig)=0.0
[253]	1530	ENDDO
	1531
	1532	c========================================================
	1533
	1534	! DO ig=1,ngridmx
	1535	! IF(tsurf(ig) .LT. 273.16-1.8) THEN
	1536	! hice(ig)=(273.16-1.8-tsurf(ig))/(273.16-1.8-240)*1.
	1537	! hice(ig)=min(hice(ig),1.0)
	1538	! ENDIF
	1539	! ENDDO
	1540
[1297]	1541
	1542
	1543
[253]	1544	c=======================================================================
[135]	1545	c Correct pressure on the new grid (menu 0)
	1546	c=======================================================================
	1547
[588]	1548
[135]	1549	if ((choix_1.eq.0).and.(.not.flagps0)) then
	1550	r = 1000.*8.31/mugaz
	1551
	1552	do j=1,jjp1
	1553	do i=1,iip1
[588]	1554	ps(i,j) = ps(i,j) *
	1555	. exp((phisold_newgrid(i,j)-phis(i,j)) /
	1556	. (t(i,j,1) * r))
[253]	1557	end do
	1558	end do
	1559
[588]	1560	c periodicite de ps en longitude
[253]	1561	do j=1,jjp1
[588]	1562	ps(1,j) = ps(iip1,j)
[135]	1563	end do
	1564	end if
	1565
[588]	1566
[135]	1567	c=======================================================================
	1568	c=======================================================================
	1569
	1570	c=======================================================================
	1571	c Initialisation de la physique / ecriture de newstartfi :
	1572	c=======================================================================
	1573
[588]	1574
[135]	1575	CALL inifilr
	1576	CALL pression(ip1jmp1, ap, bp, ps, p3d)
	1577
	1578	c-----------------------------------------------------------------------
	1579	c Initialisation pks:
	1580	c-----------------------------------------------------------------------
	1581
	1582	CALL exner_hyb(ip1jmp1, ps, p3d, beta, pks, pk, pkf)
	1583	! Calcul de la temperature potentielle teta
	1584
	1585	if (flagtset) then
	1586	DO l=1,llm
	1587	DO j=1,jjp1
	1588	DO i=1,iim
	1589	teta(i,j,l) = Tset(i,j,l) * cpp/pk(i,j,l)
	1590	ENDDO
	1591	teta (iip1,j,l)= teta (1,j,l)
	1592	ENDDO
	1593	ENDDO
	1594	else if (choix_1.eq.0) then
	1595	DO l=1,llm
	1596	DO j=1,jjp1
	1597	DO i=1,iim
	1598	teta(i,j,l) = t(i,j,l) * cpp/pk(i,j,l)
	1599	ENDDO
	1600	teta (iip1,j,l)= teta (1,j,l)
	1601	ENDDO
	1602	ENDDO
	1603	endif
	1604
	1605	C Calcul intermediaire
	1606	c
	1607	if (choix_1.eq.0) then
	1608	CALL massdair( p3d, masse )
	1609	c
	1610	print *,' ALPHAX ',alphax
	1611
	1612	DO l = 1, llm
	1613	DO i = 1, iim
	1614	xppn(i) = aire( i, 1 ) * masse( i , 1 , l )
	1615	xpps(i) = aire( i,jjp1 ) * masse( i , jjp1 , l )
	1616	ENDDO
	1617	xpn = SUM(xppn)/apoln
	1618	xps = SUM(xpps)/apols
	1619	DO i = 1, iip1
	1620	masse( i , 1 , l ) = xpn
	1621	masse( i , jjp1 , l ) = xps
	1622	ENDDO
	1623	ENDDO
	1624	endif
	1625	phis(iip1,:) = phis(1,:)
	1626
	1627	itau=0
	1628	if (choix_1.eq.0) then
	1629	day_ini=int(date)
	1630	endif
	1631	c
	1632	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
	1633
	1634	CALL caldyn0( itau,ucov,vcov,teta,ps,masse,pk,phis ,
	1635	* phi,w, pbaru,pbarv,day_ini+time )
	1636
[588]	1637
[1415]	1638	CALL dynredem0("restart.nc",day_ini,phis)
	1639	CALL dynredem1("restart.nc",0.0,vcov,ucov,teta,q,masse,ps)
[135]	1640	C
	1641	C Ecriture etat initial physique
	1642	C
	1643
[1216]	1644	call physdem0("restartfi.nc",lonfi,latfi,nsoilmx,ngridmx,llm,
	1645	& nqtot,dtphys,real(day_ini),0.0,
	1646	& airefi,albfi,ithfi,zmea,zstd,zsig,zgam,zthe)
[1415]	1647	call physdem1("restartfi.nc",nsoilmx,ngridmx,llm,nqtot,
[1216]	1648	& dtphys,real(day_ini),
	1649	& tsurf,tsoil,emis,q2,qsurf,
[1297]	1650	& cloudfrac,totalfrac,hice,
	1651	& rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
[253]	1652
[135]	1653	c=======================================================================
[1543]	1654	c Formats
[135]	1655	c=======================================================================
	1656
	1657	1 FORMAT(//10x,'la valeur de im =',i4,2x,'lue sur le fichier de dema
	1658	*rrage est differente de la valeur parametree iim =',i4//)
	1659	2 FORMAT(//10x,'la valeur de jm =',i4,2x,'lue sur le fichier de dema
	1660	*rrage est differente de la valeur parametree jjm =',i4//)
	1661	3 FORMAT(//10x,'la valeur de lllm =',i4,2x,'lue sur le fichier demar
	1662	*rage est differente de la valeur parametree llm =',i4//)
	1663
[1216]	1664	write(,) "newstart: All is well that ends well."
	1665
[135]	1666	end
	1667

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

Download in other formats:

Original Format