Context Navigation

newstart.F @ 1538

Last change on this file since 1538 was 1538, checked in by emillour, 9 years ago

Generic GCM:

Made nsoilmx be no longer a "parameter" and thus added the possibility to define the number of subsurface layers nsoilmx, along with first layer thickness "lay1_soil" and companion coefficient "alpha_soil", in callphys.def at run time. As before (when these were hard-coded), these are such that the depth of soil mid-layers are: mlayer(k)=lay1_soil*alpha_soil(k-1/2), for k=0,nsoil-1

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

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

Download in other formats:

Original Format