Context Navigation

newstart.F @ 1890

Last change on this file since 1890 was 1890, checked in by jvatant, 7 years ago
Making chemistry handling more flexible - step 2.5 + For more convenience I introduce specific modules for chemistry stuff specific to start2archive and newstart and not to pollute main module comchem_h. --JVO
File size: 40.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
[1871]	9	c Objet: Create or modify the initial state for the LMD Titan GCM
[135]	10	c ----- (fichiers NetCDF start et startfi)
	11	c
	12	c
	13	c=======================================================================
	14
[1543]	15	use mod_phys_lmdz_para, only: is_parallel, is_sequential,
	16	& is_mpi_root, is_omp_root,
	17	& is_master
[1415]	18	use infotrac, only: infotrac_init, nqtot, tname
[1890]	19	USE comchem_h, ONLY: nlaykim_up
	20	USE comchem_newstart_h, ONLY: nlaykimold
[1216]	21	USE comsoil_h, ONLY: nsoilmx, layer, mlayer, inertiedat
	22	USE surfdat_h, ONLY: phisfi, albedodat,
	23	& zmea, zstd, zsig, zgam, zthe
[1470]	24	use datafile_mod, only: datadir, surfdir
[1524]	25	use ioipsl_getin_p_mod, only: getin_p
[1415]	26	use control_mod, only: day_step, iphysiq, anneeref, planet_type
[1216]	27	use phyredem, only: physdem0, physdem1
	28	use iostart, only: open_startphy
[1403]	29	use filtreg_mod, only: inifilr
[1543]	30	USE mod_const_mpi, ONLY: COMM_LMDZ
[1422]	31	USE comvert_mod, ONLY: ap,bp,aps,bps,pa,preff
	32	USE comconst_mod, ONLY: lllm,daysec,dtvr,dtphys,cpp,kappa,
[1543]	33	. rad,omeg,g,r,pi
[1422]	34	USE serre_mod, ONLY: alphax
	35	USE temps_mod, ONLY: day_ini
	36	USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
[1670]	37	use tabfi_mod, only: tabfi
[1890]	38	use callkeys_mod, only: callchim
[1543]	39	use iniphysiq_mod, only: iniphysiq
[1670]	40	use phyetat0_mod, only: phyetat0
[1815]	41	use tracer_h
[135]	42	implicit none
	43
[1543]	44	include "dimensions.h"
[1308]	45	integer, parameter :: ngridmx = (2+(jjm-1)*iim - 1/jjm)
[1543]	46	include "paramet.h"
	47	include "comgeom2.h"
	48	include "comdissnew.h"
	49	include "netcdf.inc"
	50
[135]	51	c=======================================================================
	52	c Declarations
	53	c=======================================================================
	54
	55	c Variables dimension du fichier "start_archive"
	56	c------------------------------------
[1543]	57	CHARACTER relief*3
[135]	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------------------
[1543]	97	REAL tsurf(ngridmx) ! surface temperature
	98	REAL,ALLOCATABLE :: tsoil(:,:) ! soil temperature
	99	REAL emis(ngridmx) ! surface emissivity
[135]	100	real emisread ! added by RW
[1216]	101	REAL,ALLOCATABLE :: qsurf(:,:)
[1308]	102	REAL q2(ngridmx,llm+1)
[135]	103	! REAL rnaturfi(ngridmx)
	104	real alb(iip1,jjp1),albfi(ngridmx) ! albedos
[1543]	105	real,ALLOCATABLE :: ith(:,:,:),ithfi(:,:) ! thermal inertia (3D)
[135]	106	real surfith(iip1,jjp1),surfithfi(ngridmx) ! surface thermal inertia (2D)
	107	REAL latfi(ngridmx),lonfi(ngridmx),airefi(ngridmx)
	108
	109	INTEGER i,j,l,isoil,ig,idum
	110	real mugaz ! molar mass of the atmosphere
	111
[1647]	112	integer ierr
[135]	113
	114	c Variables on the new grid along scalar points
	115	c------------------------------------------------------
	116	! REAL p(iip1,jjp1)
	117	REAL t(iip1,jjp1,llm)
	118	REAL tset(iip1,jjp1,llm)
	119	real phisold_newgrid(iip1,jjp1)
	120	REAL :: teta(iip1, jjp1, llm)
	121	REAL :: pk(iip1,jjp1,llm)
	122	REAL :: pkf(iip1,jjp1,llm)
	123	REAL :: ps(iip1, jjp1)
	124	REAL :: masse(iip1,jjp1,llm)
	125	REAL :: xpn,xps,xppn(iim),xpps(iim)
	126	REAL :: p3d(iip1, jjp1, llm+1)
	127	REAL :: beta(iip1,jjp1,llm)
	128	! REAL dteta(ip1jmp1,llm)
	129
	130	c Variable de l'ancienne grille
	131	c------------------------------
	132	real time
	133	real tab_cntrl(100)
	134	real tab_cntrl_bis(100)
[1889]	135
[135]	136	c variables diverses
	137	c-------------------
	138	real choix_1,pp
	139	character*80 fichnom
[988]	140	character*250 filestring
[1543]	141	integer Lmodif,iq
[135]	142	character modif*20
	143	real z_reel(iip1,jjp1)
	144	real tsud,albsud,alb_bb,ith_bb,Tiso,Tabove
	145	real ptoto,pcap,patm,airetot,ptotn,patmn,psea
	146	! real ssum
	147	character*1 yes
	148	logical :: flagtset=.false. , flagps0=.false.
[1369]	149	real val, val2, val3, val4 ! to store temporary variables
[135]	150
	151	INTEGER :: itau
	152
	153	character(len=20) :: txt ! to store some text
[988]	154	character(len=50) :: surfacefile ! "surface.nc" (or equivalent file)
	155	character(len=150) :: longline
[135]	156	integer :: count
[699]	157	real :: profile(llm+1) ! to store an atmospheric profile + surface value
[135]	158
[253]	159	! added by BC for equilibrium temperature startup
	160	real teque
	161
	162	! added by RW for nuketharsis
	163	real fact1
	164	real fact2
	165
[1789]	166	! added by JVO for Titan
	167	REAL tankCH4(ngridmx)
[253]	168
[135]	169	c sortie visu pour les champs dynamiques
	170	c---------------------------------------
	171	! INTEGER :: visuid
	172	! real :: time_step,t_ops,t_wrt
	173	! CHARACTER*80 :: visu_file
	174
[535]	175	cpp = 0.
	176	preff = 0.
	177	pa = 0. ! to ensure disaster rather than minor error if we don`t
	178	! make deliberate choice of these values elsewhere.
[135]	179
[1644]	180	planet_type="titan"
[1543]	181
[1216]	182	! initialize "serial/parallel" related stuff
[1543]	183	! (required because we call tabfi() below, before calling iniphysiq)
	184	is_sequential=.true.
	185	is_parallel=.false.
	186	is_mpi_root=.true.
	187	is_omp_root=.true.
	188	is_master=.true.
	189
[1216]	190	! Load tracer number and names:
[1415]	191	call infotrac_init
[1216]	192	! allocate arrays
	193	allocate(q(iip1,jjp1,llm,nqtot))
	194	allocate(qsurf(ngridmx,nqtot))
[1415]	195
[1543]	196	! get value of nsoilmx and allocate corresponding arrays
[1589]	197	! a default value of nsoilmx is set in comsoil_h
[1543]	198	call getin_p("nsoilmx",nsoilmx)
	199
	200	allocate(tsoil(ngridmx,nsoilmx))
	201	allocate(ith(iip1,jjp1,nsoilmx),ithfi(ngridmx,nsoilmx))
	202
[135]	203	c=======================================================================
	204	c Choice of the start file(s) to use
	205	c=======================================================================
	206	write(,) 'From which kind of files do you want to create new',
	207	. 'start and startfi files'
	208	write(,) ' 0 - from a file start_archive'
	209	write(,) ' 1 - from files start and startfi'
	210
	211	c-----------------------------------------------------------------------
	212	c Open file(s) to modify (start or start_archive)
	213	c-----------------------------------------------------------------------
	214
	215	DO
	216	read(,,iostat=ierr) choix_1
	217	if ((choix_1 /= 0).OR.(choix_1 /=1)) EXIT
	218	ENDDO
	219
	220	c Open start_archive
	221	c ~~~~~~~~~~~~~~~~~~~~~~~~~~
	222	if (choix_1.eq.0) then
	223
	224	write(,) 'Creating start files from:'
	225	write(,) './start_archive.nc'
	226	write(,)
	227	fichnom = 'start_archive.nc'
	228	ierr = NF_OPEN (fichnom, NF_NOWRITE,nid)
	229	IF (ierr.NE.NF_NOERR) THEN
	230	write(6,*)' Problem opening file:',fichnom
	231	write(6,*)' ierr = ', ierr
	232	CALL ABORT
	233	ENDIF
	234	tab0 = 50
	235	Lmodif = 1
	236
	237	c OR open start and startfi files
	238	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	239	else
	240	write(,) 'Creating start files from:'
	241	write(,) './start.nc and ./startfi.nc'
	242	write(,)
	243	fichnom = 'start.nc'
	244	ierr = NF_OPEN (fichnom, NF_NOWRITE,nid_dyn)
	245	IF (ierr.NE.NF_NOERR) THEN
	246	write(6,*)' Problem opening file:',fichnom
	247	write(6,*)' ierr = ', ierr
	248	CALL ABORT
	249	ENDIF
	250
	251	fichnom = 'startfi.nc'
	252	ierr = NF_OPEN (fichnom, NF_NOWRITE,nid_fi)
	253	IF (ierr.NE.NF_NOERR) THEN
	254	write(6,*)' Problem opening file:',fichnom
	255	write(6,*)' ierr = ', ierr
	256	CALL ABORT
	257	ENDIF
	258
	259	tab0 = 0
	260	Lmodif = 0
	261
	262	endif
	263
[649]	264	c=======================================================================
	265	c INITIALISATIONS DIVERSES
	266	c=======================================================================
	267
[787]	268	! Initialize global tracer indexes (stored in tracer.h)
	269	! ... this has to be done before phyetat0
[1815]	270	call initracer2(nqtot,tname)
[787]	271
[837]	272	! Take care of arrays in common modules
	273	! ALLOCATE ARRAYS in surfdat_h (if not already done, e.g. when using start_archive)
	274	IF (.not. ALLOCATED(albedodat)) ALLOCATE(albedodat(ngridmx))
	275	IF (.not. ALLOCATED(phisfi)) ALLOCATE(phisfi(ngridmx))
	276	IF (.not. ALLOCATED(zmea)) ALLOCATE(zmea(ngridmx))
	277	IF (.not. ALLOCATED(zstd)) ALLOCATE(zstd(ngridmx))
	278	IF (.not. ALLOCATED(zsig)) ALLOCATE(zsig(ngridmx))
	279	IF (.not. ALLOCATED(zgam)) ALLOCATE(zgam(ngridmx))
	280	IF (.not. ALLOCATED(zthe)) ALLOCATE(zthe(ngridmx))
[787]	281
[135]	282	c-----------------------------------------------------------------------
	283	c Lecture du tableau des parametres du run (pour la dynamique)
	284	c-----------------------------------------------------------------------
	285	if (choix_1.eq.0) then
	286
	287	write(,) 'reading tab_cntrl START_ARCHIVE'
	288	c
	289	ierr = NF_INQ_VARID (nid, "controle", nvarid)
	290	#ifdef NC_DOUBLE
	291	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tab_cntrl)
	292	#else
	293	ierr = NF_GET_VAR_REAL(nid, nvarid, tab_cntrl)
	294	#endif
	295	c
	296	else if (choix_1.eq.1) then
	297
	298	write(,) 'reading tab_cntrl START'
	299	c
	300	ierr = NF_INQ_VARID (nid_dyn, "controle", nvarid)
	301	#ifdef NC_DOUBLE
	302	ierr = NF_GET_VAR_DOUBLE(nid_dyn, nvarid, tab_cntrl)
	303	#else
	304	ierr = NF_GET_VAR_REAL(nid_dyn, nvarid, tab_cntrl)
	305	#endif
	306	c
	307	write(,) 'reading tab_cntrl STARTFI'
	308	c
	309	ierr = NF_INQ_VARID (nid_fi, "controle", nvarid)
	310	#ifdef NC_DOUBLE
	311	ierr = NF_GET_VAR_DOUBLE(nid_fi, nvarid, tab_cntrl_bis)
	312	#else
	313	ierr = NF_GET_VAR_REAL(nid_fi, nvarid, tab_cntrl_bis)
	314	#endif
	315	c
	316	do i=1,50
	317	tab_cntrl(i+50)=tab_cntrl_bis(i)
	318	enddo
[649]	319	write(,) 'printing tab_cntrl', tab_cntrl
	320	do i=1,100
	321	write(,) i,tab_cntrl(i)
	322	enddo
[1543]	323
[649]	324	write(,) 'Reading file START'
	325	fichnom = 'start.nc'
[1421]	326	CALL dynetat0(fichnom,vcov,ucov,teta,q,masse,
[649]	327	. ps,phis,time)
	328
[1589]	329	CALL iniconst
	330	CALL inigeom
	331
	332	! Initialize the physics
	333	CALL iniphysiq(iim,jjm,llm,
	334	& (jjm-1)*iim+2,comm_lmdz,
	335	& daysec,day_ini,dtphys,
	336	& rlatu,rlatv,rlonu,rlonv,
	337	& aire,cu,cv,rad,g,r,cpp,
	338	& 1)
	339
	340	! Lmodif set to 0 to disable modifications possibility in phyeta0
	341	Lmodif=0
[649]	342	write(,) 'Reading file STARTFI'
	343	fichnom = 'startfi.nc'
[1670]	344	CALL phyetat0(.true.,ngridmx,llm,fichnom,tab0,Lmodif,nsoilmx,
[1308]	345	. nqtot,day_ini,time,
[1789]	346	. tsurf,tsoil,emis,q2,qsurf,tankCH4)
[649]	347
[905]	348	! copy albedo and soil thermal inertia on (local) physics grid
[649]	349	do i=1,ngridmx
	350	albfi(i) = albedodat(i)
[1543]	351	do j=1,nsoilmx
[649]	352	ithfi(i,j) = inertiedat(i,j)
[1543]	353	enddo
[649]	354	! build a surfithfi(:) using 1st layer of ithfi(:), which might
[905]	355	! be needed later on if reinitializing soil thermal inertia
[649]	356	surfithfi(i)=ithfi(i,1)
	357	enddo
[905]	358	! also copy albedo and soil thermal inertia on (local) dynamics grid
	359	! so that options below can manipulate either (but must then ensure
	360	! to correctly recast things on physics grid)
	361	call gr_fi_dyn(1,ngridmx,iip1,jjp1,albfi,alb)
	362	call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
	363	call gr_fi_dyn(1,ngridmx,iip1,jjp1,surfithfi,surfith)
[135]	364
	365	endif
	366	c-----------------------------------------------------------------------
[1543]	367	c Initialisation des constantes dynamique
[135]	368	c-----------------------------------------------------------------------
	369
[1850]	370	! JVO 2017 : Added a +1 shift in tab_cntrl indexes to be coherent
	371	! with Earth-like Titan start files ( cf dynetat0 from common )
[1849]	372
	373	kappa = tab_cntrl(10)
	374	etot0 = tab_cntrl(13)
	375	ptot0 = tab_cntrl(14)
	376	ztot0 = tab_cntrl(15)
	377	stot0 = tab_cntrl(16)
	378	ang0 = tab_cntrl(17)
[135]	379	write(,) "Newstart: kappa,etot0,ptot0,ztot0,stot0,ang0"
	380	write(,) kappa,etot0,ptot0,ztot0,stot0,ang0
	381
	382	! for vertical coordinate
[1849]	383	preff=tab_cntrl(19) ! reference surface pressure
	384	pa=tab_cntrl(18) ! reference pressure at which coord is purely pressure
[135]	385	!NB: in start_archive files tab_cntrl(17)=tab_cntrl(18)=0
	386	write(,) "Newstart: preff=",preff," pa=",pa
	387	yes=' '
	388	do while ((yes.ne.'y').and.(yes.ne.'n'))
	389	write(,) "Change the values of preff and pa? (y/n)"
	390	read(*,fmt='(a)') yes
	391	end do
	392	if (yes.eq.'y') then
	393	write(,)"New value of reference surface pressure preff?"
[1849]	394	write(,)" (for Titan, typically preff=140000)"
[135]	395	read(,) preff
	396	write(,)"New value of reference pressure pa for purely"
	397	write(,)"pressure levels (for hybrid coordinates)?"
[1849]	398	write(,)" (for Titan, typically pa=10000)"
[135]	399	read(,) pa
	400	endif
	401	c-----------------------------------------------------------------------
	402	c Lecture du tab_cntrl et initialisation des constantes physiques
	403	c - pour start: Lmodif = 0 => pas de modifications possibles
	404	c (modif dans le tabfi de readfi + loin)
	405	c - pour start_archive: Lmodif = 1 => modifications possibles
	406	c-----------------------------------------------------------------------
	407	if (choix_1.eq.0) then
[1216]	408	! tabfi requires that input file be first opened by open_startphy(fichnom)
	409	fichnom = 'start_archive.nc'
	410	call open_startphy(fichnom)
[787]	411	call tabfi (ngridmx,nid,Lmodif,tab0,day_ini,lllm,p_rad,
[135]	412	. p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
	413	else if (choix_1.eq.1) then
[1216]	414	fichnom = 'startfi.nc'
	415	call open_startphy(fichnom)
[649]	416	Lmodif=1 ! Lmodif set to 1 to allow modifications in phyeta0
[787]	417	call tabfi (ngridmx,nid_fi,Lmodif,tab0,day_ini,lllm,p_rad,
[135]	418	. p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
	419	endif
	420
[1321]	421	if (p_omeg .eq. -9999.) then
	422	p_omeg = 8.*atan(1.)/p_daysec
	423	print*,"new value of omega",p_omeg
	424	endif
	425
[135]	426	rad = p_rad
	427	omeg = p_omeg
	428	g = p_g
	429	cpp = p_cpp
	430	mugaz = p_mugaz
	431	daysec = p_daysec
	432
[1330]	433	if (p_omeg .eq. -9999.) then
	434	p_omeg = 8.*atan(1.)/p_daysec
	435	print*,"new value of omega",p_omeg
	436	endif
	437
[535]	438	kappa = 8.3141000./(p_mugazp_cpp) ! added by RDW
[135]	439
	440	c=======================================================================
	441	c INITIALISATIONS DIVERSES
	442	c=======================================================================
	443	! Load parameters from run.def file
	444	CALL defrun_new( 99, .TRUE. )
[1589]	445	! Initialize dynamics geometry and co. (which, when using start.nc may
	446	! have changed because of modifications (tabi, preff, pa) above)
[135]	447	CALL iniconst
	448	CALL inigeom
	449	idum=-1
	450	idum=0
	451
[1589]	452	! Initialize the physics for start_archive only
	453	IF (choix_1.eq.0) THEN
[1543]	454	CALL iniphysiq(iim,jjm,llm,
	455	& (jjm-1)*iim+2,comm_lmdz,
	456	& daysec,day_ini,dtphys,
	457	& rlatu,rlatv,rlonu,rlonv,
	458	& aire,cu,cv,rad,g,r,cpp,
	459	& 1)
[1589]	460	ENDIF
[135]	461
	462	c=======================================================================
	463	c lecture topographie, albedo, inertie thermique, relief sous-maille
	464	c=======================================================================
	465
[988]	466	if (choix_1.eq.0) then ! for start_archive files,
	467	! where an external "surface.nc" file is needed
[135]	468
	469	c do while((relief(1:3).ne.'mol').AND.(relief(1:3).ne.'pla'))
	470	c write(,)
	471	c write(,) 'choix du relief (mola,pla)'
	472	c write(,) '(Topographie MGS MOLA, plat)'
	473	c read(*,fmt='(a3)') relief
	474	relief="mola"
	475	c enddo
	476
[988]	477	! First get the correct value of datadir, if not already done:
	478	! default 'datadir' is set in "datafile_mod"
[1315]	479	call getin_p("datadir",datadir)
[988]	480	write(,) 'Available surface data files are:'
[1470]	481	filestring='ls '//trim(datadir)//'/'//
	482	& trim(surfdir)//' \| grep .nc'
[988]	483	call system(filestring)
[1470]	484	! but in ye old days these files were in datadir, so scan it as well
	485	! for the sake of retro-compatibility
	486	filestring='ls '//trim(datadir)//'/'//' \| grep .nc'
	487	call system(filestring)
[135]	488
[988]	489	write(,) ''
	490	write(,) 'Please choose the relevant file',
	491	& ' (e.g. "surface_mars.nc")'
	492	write(,) ' or "none" to not use any (and set planetary'
	493	write(,) ' albedo and surface thermal inertia)'
	494	read(*,fmt='(a50)') surfacefile
[135]	495
[988]	496	if (surfacefile.ne."none") then
[135]	497
[988]	498	CALL datareadnc(relief,surfacefile,phis,alb,surfith,
	499	& zmeaS,zstdS,zsigS,zgamS,ztheS)
	500	else
	501	! specific case when not using a "surface.nc" file
	502	phis(:,:)=0
	503	zmeaS(:,:)=0
	504	zstdS(:,:)=0
	505	zsigS(:,:)=0
	506	zgamS(:,:)=0
	507	ztheS(:,:)=0
	508
	509	write(,) "Enter value of albedo of the bare ground:"
	510	read(,) alb(1,1)
	511	alb(:,:)=alb(1,1)
	512
	513	write(,) "Enter value of thermal inertia of soil:"
	514	read(,) surfith(1,1)
	515	surfith(:,:)=surfith(1,1)
	516
	517	endif ! of if (surfacefile.ne."none")
[135]	518
[988]	519	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
	520	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,surfith,surfithfi)
	521	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
	522	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zmeaS,zmea)
	523	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zstdS,zstd)
	524	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zsigS,zsig)
	525	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zgamS,zgam)
	526	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ztheS,zthe)
	527
	528	endif ! of if (choix_1.eq.0)
	529
	530
[135]	531	c=======================================================================
	532	c Lecture des fichiers (start ou start_archive)
	533	c=======================================================================
	534
	535	if (choix_1.eq.0) then
	536
	537	write(,) 'Reading file START_ARCHIVE'
[1308]	538	CALL lect_start_archive(ngridmx,llm,
	539	& date,tsurf,tsoil,emis,q2,
[1889]	540	& t,ucov,vcov,ps,teta,phisold_newgrid,
	541	& q,qsurf,surfith,nid)
[135]	542	write(,) "OK, read start_archive file"
[1543]	543	! copy soil thermal inertia
	544	ithfi(:,:)=inertiedat(:,:)
	545
[135]	546	ierr= NF_CLOSE(nid)
	547
[649]	548	else if (choix_1.eq.1) then
	549	!do nothing, start and startfi have already been read
[135]	550	else
	551	CALL exit(1)
	552	endif
	553
	554	dtvr = daysec/FLOAT(day_step)
	555	dtphys = dtvr * FLOAT(iphysiq)
	556
	557	c=======================================================================
	558	c
	559	c=======================================================================
	560
	561	do ! infinite loop on list of changes
	562
	563	write(,)
	564	write(,)
	565	write(,) 'List of possible changes :'
	566	write(,) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
	567	write(,)
	568	write(,) 'flat : no topography ("aquaplanet")'
[1359]	569	write(,) 'set_ps_to_preff : used if changing preff with topo'
[135]	570	write(,) 'bilball : uniform albedo and thermal inertia'
	571	write(,) 'qname : change tracer name'
[996]	572	write(,) 't=profile : read temperature profile in profile.in'
[135]	573	write(,) 'q=0 : ALL tracer =zero'
	574	write(,) 'q=x : give a specific uniform value to one tracer'
[699]	575	write(,) 'q=profile : specify a profile for a tracer'
[253]	576	write(,) 'isotherm : Isothermal Temperatures, wind set to zero'
	577	write(,) 'radequi : Earth-like radiative equilibrium temperature
	578	$ profile (lat-alt) and winds set to zero'
[135]	579	write(,) 'ptot : change total pressure'
	580	write(,) 'emis : change surface emissivity'
[253]	581	write(,) 'therm_ini_s : Set soil thermal inertia to reference sur
[135]	582	&face values'
[1297]	583	write(,) 'slab_ocean_0 : initialisation of slab
	584	$ocean variables'
[135]	585
	586	write(,)
	587	write(,) 'Change to perform ?'
	588	write(,) ' (enter keyword or return to end)'
	589	write(,)
	590
	591	read(*,fmt='(a20)') modif
[1216]	592	if (modif(1:1) .eq. ' ')then
	593	exit ! exit loop on changes
	594	endif
[135]	595
	596	write(,)
	597	write(,) trim(modif) , ' : '
	598
	599	c 'flat : no topography ("aquaplanet")'
	600	c -------------------------------------
[699]	601	if (trim(modif) .eq. 'flat') then
[135]	602	c set topo to zero
[699]	603	z_reel(1:iip1,1:jjp1)=0
	604	phis(1:iip1,1:jjp1)=g*z_reel(1:iip1,1:jjp1)
[135]	605	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
	606	write(,) 'topography set to zero.'
	607	write(,) 'WARNING : the subgrid topography parameters',
	608	& ' were not set to zero ! => set calllott to F'
	609
[253]	610	c Choice of surface pressure
[135]	611	yes=' '
	612	do while ((yes.ne.'y').and.(yes.ne.'n'))
	613	write(,) 'Do you wish to choose homogeneous surface',
	614	& 'pressure (y) or let newstart interpolate ',
	615	& ' the previous field (n)?'
	616	read(*,fmt='(a)') yes
	617	end do
	618	if (yes.eq.'y') then
	619	flagps0=.true.
	620	write(,) 'New value for ps (Pa) ?'
	621	201 read(,,iostat=ierr) patm
	622	if(ierr.ne.0) goto 201
[1320]	623	write(,) patm
	624	if (patm.eq.-9999.) then
	625	patm = preff
	626	write(,) "we set patm = preff", preff
	627	endif
[135]	628	write(,)
	629	write(,) ' new ps everywhere (Pa) = ', patm
	630	write(,)
	631	do j=1,jjp1
	632	do i=1,iip1
	633	ps(i,j)=patm
	634	enddo
	635	enddo
	636	end if
	637
[1359]	638	c 'set_ps_to_preff' : used if changing preff with topo
	639	c ----------------------------------------------------
	640	else if (trim(modif) .eq. 'set_ps_to_preff') then
	641	do j=1,jjp1
	642	do i=1,iip1
	643	ps(i,j)=preff
	644	enddo
	645	enddo
	646
[253]	647
	648	c bilball : uniform albedo, thermal inertia
	649	c -----------------------------------------
[699]	650	else if (trim(modif) .eq. 'bilball') then
[135]	651	write(,) 'constante albedo and iner.therm:'
	652	write(,) 'New value for albedo (ex: 0.25) ?'
	653	101 read(,,iostat=ierr) alb_bb
	654	if(ierr.ne.0) goto 101
	655	write(,)
	656	write(,) ' uniform albedo (new value):',alb_bb
	657	write(,)
	658
	659	write(,) 'New value for thermal inertia (eg: 247) ?'
	660	102 read(,,iostat=ierr) ith_bb
	661	if(ierr.ne.0) goto 102
	662	write(,) 'uniform thermal inertia (new value):',ith_bb
	663	DO j=1,jjp1
	664	DO i=1,iip1
[1543]	665	alb(i,j) = alb_bb ! albedo
	666	do isoil=1,nsoilmx
	667	ith(i,j,isoil) = ith_bb ! thermal inertia
	668	enddo
[135]	669	END DO
	670	END DO
	671	! CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ith,ithfi)
	672	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
	673	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
	674
	675	c ptot : Modification of the total pressure: ice + current atmosphere
	676	c -------------------------------------------------------------------
[699]	677	else if (trim(modif).eq.'ptot') then
[135]	678
	679	c calcul de la pression totale glace + atm actuelle
	680	patm=0.
	681	airetot=0.
	682	pcap=0.
	683	DO j=1,jjp1
	684	DO i=1,iim
	685	ig=1+(j-2)*iim +i
	686	if(j.eq.1) ig=1
	687	if(j.eq.jjp1) ig=ngridmx
	688	patm = patm + ps(i,j)*aire(i,j)
	689	airetot= airetot + aire(i,j)
	690	ENDDO
	691	ENDDO
	692	ptoto = pcap + patm
	693
[1647]	694	print*,'Current total pressure at surface (atm) ',
[135]	695	& ptoto/airetot
	696
	697	print*,'new value?'
	698	read(,) ptotn
	699	ptotn=ptotn*airetot
	700	patmn=ptotn-pcap
	701	print*,'ptoto,patm,ptotn,patmn'
	702	print*,ptoto,patm,ptotn,patmn
	703	print*,'Mult. factor for pressure (atm only)', patmn/patm
	704	do j=1,jjp1
	705	do i=1,iip1
	706	ps(i,j)=ps(i,j)*patmn/patm
	707	enddo
	708	enddo
	709
[253]	710
	711
[135]	712	c Correction pour la conservation des traceurs
	713	yes=' '
	714	do while ((yes.ne.'y').and.(yes.ne.'n'))
	715	write(,) 'Do you wish to conserve tracer total mass (y)',
	716	& ' or tracer mixing ratio (n) ?'
	717	read(*,fmt='(a)') yes
	718	end do
	719
	720	if (yes.eq.'y') then
	721	write(,) 'OK : conservation of tracer total mass'
[1216]	722	DO iq =1, nqtot
[135]	723	DO l=1,llm
	724	DO j=1,jjp1
	725	DO i=1,iip1
	726	q(i,j,l,iq)=q(i,j,l,iq)*patm/patmn
	727	ENDDO
	728	ENDDO
	729	ENDDO
	730	ENDDO
	731	else
	732	write(,) 'OK : conservation of tracer mixing ratio'
	733	end if
	734
	735	c Correction pour la pression au niveau de la mer
	736	yes=' '
	737	do while ((yes.ne.'y').and.(yes.ne.'n'))
	738	write(,) 'Do you wish fix pressure at sea level (y)',
	739	& ' or not (n) ?'
	740	read(*,fmt='(a)') yes
	741	end do
	742
	743	if (yes.eq.'y') then
	744	write(,) 'Value?'
[1543]	745	read(,,iostat=ierr) psea
[135]	746	DO i=1,iip1
	747	DO j=1,jjp1
	748	ps(i,j)=psea
	749
	750	ENDDO
	751	ENDDO
[1543]	752	write(,) 'psea=',psea
[135]	753	else
	754	write(,) 'no'
	755	end if
	756
	757
	758	c emis : change surface emissivity (added by RW)
	759	c ----------------------------------------------
	760	else if (trim(modif).eq.'emis') then
	761
	762	print*,'new value?'
	763	read(,) emisread
	764
	765	do i=1,ngridmx
	766	emis(i)=emisread
	767	enddo
	768
	769	c qname : change tracer name
	770	c --------------------------
	771	else if (trim(modif).eq.'qname') then
	772	yes='y'
	773	do while (yes.eq.'y')
	774	write(,) 'Which tracer name do you want to change ?'
[1216]	775	do iq=1,nqtot
	776	write(*,'(i3,a3,a20)')iq,' : ',trim(tname(iq))
[135]	777	enddo
	778	write(*,'(a35,i3)')
[1216]	779	& '(enter tracer number; between 1 and ',nqtot
[135]	780	write(,)' or any other value to quit this option)'
	781	read(,) iq
[1216]	782	if ((iq.ge.1).and.(iq.le.nqtot)) then
	783	write(,)'Change tracer name ',trim(tname(iq)),' to ?'
[135]	784	read(,) txt
[1216]	785	tname(iq)=txt
[135]	786	write(,)'Do you want to change another tracer name (y/n)?'
	787	read(*,'(a)') yes
	788	else
	789	! inapropiate value of iq; quit this option
	790	yes='n'
[1216]	791	endif ! of if ((iq.ge.1).and.(iq.le.nqtot))
[135]	792	enddo ! of do while (yes.ne.'y')
	793
	794	c q=0 : set tracers to zero
	795	c -------------------------
[699]	796	else if (trim(modif).eq.'q=0') then
[135]	797	c mise a 0 des q (traceurs)
	798	write(,) 'Tracers set to 0 (1.E-30 in fact)'
[1216]	799	DO iq =1, nqtot
[135]	800	DO l=1,llm
	801	DO j=1,jjp1
	802	DO i=1,iip1
	803	q(i,j,l,iq)=1.e-30
	804	ENDDO
	805	ENDDO
	806	ENDDO
	807	ENDDO
	808
	809	c set surface tracers to zero
[1216]	810	DO iq =1, nqtot
[135]	811	DO ig=1,ngridmx
	812	qsurf(ig,iq)=0.
	813	ENDDO
	814	ENDDO
	815
	816	c q=x : initialise tracer manually
	817	c --------------------------------
[699]	818	else if (trim(modif).eq.'q=x') then
[135]	819	write(,) 'Which tracer do you want to modify ?'
[1216]	820	do iq=1,nqtot
	821	write(,)iq,' : ',trim(tname(iq))
[135]	822	enddo
[1216]	823	write(,) '(choose between 1 and ',nqtot,')'
[135]	824	read(,) iq
[1216]	825	write(,)'mixing ratio of tracer ',trim(tname(iq)),
[135]	826	& ' ? (kg/kg)'
	827	read(,) val
	828	DO l=1,llm
	829	DO j=1,jjp1
	830	DO i=1,iip1
	831	q(i,j,l,iq)=val
	832	ENDDO
	833	ENDDO
	834	ENDDO
[1216]	835	write(,) 'SURFACE value of tracer ',trim(tname(iq)),
[135]	836	& ' ? (kg/m2)'
	837	read(,) val
	838	DO ig=1,ngridmx
	839	qsurf(ig,iq)=val
	840	ENDDO
	841
[996]	842	c t=profile : initialize temperature with a given profile
	843	c -------------------------------------------------------
	844	else if (trim(modif) .eq. 't=profile') then
	845	write(,) 'Temperature profile from ASCII file'
	846	write(,) "'profile.in' e.g. 1D output"
	847	write(,) "(one value per line in file; starting with"
	848	write(,) "surface value, the 1st atmospheric layer"
	849	write(,) "followed by 2nd, etc. up to top of atmosphere)"
	850	txt="profile.in"
	851	open(33,file=trim(txt),status='old',form='formatted',
	852	& iostat=ierr)
	853	if (ierr.eq.0) then
	854	! OK, found file 'profile_...', load the profile
	855	do l=1,llm+1
	856	read(33,*,iostat=ierr) profile(l)
	857	write(,) profile(l)
	858	if (ierr.ne.0) then ! something went wrong
	859	exit ! quit loop
	860	endif
	861	enddo
	862	if (ierr.eq.0) then
	863	tsurf(1:ngridmx)=profile(1)
	864	tsoil(1:ngridmx,1:nsoilmx)=profile(1)
	865	do l=1,llm
	866	Tset(1:iip1,1:jjp1,l)=profile(l+1)
	867	flagtset=.true.
	868	enddo
	869	ucov(1:iip1,1:jjp1,1:llm)=0.
	870	vcov(1:iip1,1:jjm,1:llm)=0.
[1308]	871	q2(1:ngridmx,1:llm+1)=0.
[996]	872	else
	873	write(,)'problem reading file ',trim(txt),' !'
	874	write(,)'No modifications to temperature'
	875	endif
	876	else
	877	write(,)'Could not find file ',trim(txt),' !'
	878	write(,)'No modifications to temperature'
	879	endif
	880
[699]	881	c q=profile : initialize tracer with a given profile
	882	c --------------------------------------------------
	883	else if (trim(modif) .eq. 'q=profile') then
	884	write(,) 'Tracer profile will be sought in ASCII file'
	885	write(,) "'profile_tracer' where 'tracer' is tracer name"
	886	write(,) "(one value per line in file; starting with"
	887	write(,) "surface value, the 1st atmospheric layer"
	888	write(,) "followed by 2nd, etc. up to top of atmosphere)"
	889	write(,) 'Which tracer do you want to set?'
[1216]	890	do iq=1,nqtot
	891	write(,)iq,' : ',trim(tname(iq))
[699]	892	enddo
[1216]	893	write(,) '(choose between 1 and ',nqtot,')'
[699]	894	read(,) iq
[1216]	895	if ((iq.lt.1).or.(iq.gt.nqtot)) then
[699]	896	! wrong value for iq, go back to menu
	897	write(,) "wrong input value:",iq
	898	cycle
	899	endif
	900	! look for input file 'profile_tracer'
[1216]	901	txt="profile_"//trim(tname(iq))
[699]	902	open(41,file=trim(txt),status='old',form='formatted',
	903	& iostat=ierr)
	904	if (ierr.eq.0) then
	905	! OK, found file 'profile_...', load the profile
	906	do l=1,llm+1
	907	read(41,*,iostat=ierr) profile(l)
	908	if (ierr.ne.0) then ! something went wrong
	909	exit ! quit loop
	910	endif
	911	enddo
	912	if (ierr.eq.0) then
	913	! initialize tracer values
	914	qsurf(:,iq)=profile(1)
	915	do l=1,llm
	916	q(:,:,l,iq)=profile(l+1)
	917	enddo
[1216]	918	write(,)'OK, tracer ',trim(tname(iq)),' initialized '
	919	& ,'using values from file ',trim(txt)
[699]	920	else
	921	write(,)'problem reading file ',trim(txt),' !'
[1216]	922	write(,)'No modifications to tracer ',trim(tname(iq))
[699]	923	endif
	924	else
	925	write(,)'Could not find file ',trim(txt),' !'
[1216]	926	write(,)'No modifications to tracer ',trim(tname(iq))
[699]	927	endif
	928
[135]	929
	930	c isotherm : Isothermal temperatures and no winds
[253]	931	c -----------------------------------------------
[699]	932	else if (trim(modif) .eq. 'isotherm') then
[135]	933
	934	write(,)'Isothermal temperature of the atmosphere,
	935	& surface and subsurface'
	936	write(,) 'Value of this temperature ? :'
	937	203 read(,,iostat=ierr) Tiso
	938	if(ierr.ne.0) goto 203
	939
[699]	940	tsurf(1:ngridmx)=Tiso
	941
	942	tsoil(1:ngridmx,1:nsoilmx)=Tiso
	943
	944	Tset(1:iip1,1:jjp1,1:llm)=Tiso
[135]	945	flagtset=.true.
[1359]	946
	947	t(1:iip1,1:jjp1,1:llm)=Tiso
	948	!! otherwise hydrost. integrations below
	949	!! use the wrong temperature
	950	!! -- NB: Tset might be useless
	951
[699]	952	ucov(1:iip1,1:jjp1,1:llm)=0
	953	vcov(1:iip1,1:jjm,1:llm)=0
[1308]	954	q2(1:ngridmx,1:llm+1)=0
[135]	955
[253]	956	c radequi : Radiative equilibrium profile of temperatures and no winds
	957	c --------------------------------------------------------------------
[699]	958	else if (trim(modif) .eq. 'radequi') then
[135]	959
[253]	960	write(,)'radiative equilibrium temperature profile'
	961
	962	do ig=1, ngridmx
	963	teque= 335.0-60.0sin(latfi(ig))sin(latfi(ig))-
	964	& 10.0cos(latfi(ig))cos(latfi(ig))
	965	tsurf(ig) = MAX(220.0,teque)
	966	end do
	967	do l=2,nsoilmx
	968	do ig=1, ngridmx
	969	tsoil(ig,l) = tsurf(ig)
	970	end do
	971	end do
	972	DO j=1,jjp1
	973	DO i=1,iim
	974	Do l=1,llm
	975	teque=335.-60.0sin(rlatu(j))sin(rlatu(j))-
	976	& 10.0cos(rlatu(j))cos(rlatu(j))
	977	Tset(i,j,l)=MAX(220.0,teque)
	978	end do
	979	end do
	980	end do
	981	flagtset=.true.
[699]	982	ucov(1:iip1,1:jjp1,1:llm)=0
	983	vcov(1:iip1,1:jjm,1:llm)=0
[1308]	984	q2(1:ngridmx,1:llm+1)=0
[253]	985
[135]	986	! therm_ini_s: (re)-set soil thermal inertia to reference surface values
	987	! ----------------------------------------------------------------------
	988
[1543]	989	else if (trim(modif) .eq. 'therm_ini_s') then
[135]	990	! write(,)"surfithfi(1):",surfithfi(1)
[1543]	991	do isoil=1,nsoilmx
	992	inertiedat(1:ngridmx,isoil)=surfithfi(1:ngridmx)
	993	enddo
[135]	994	write(,)'OK: Soil thermal inertia has been reset to referenc
	995	&e surface values'
[1543]	996	! write(,)"inertiedat(1,1):",inertiedat(1,1)
	997	ithfi(:,:)=inertiedat(:,:)
	998	! recast ithfi() onto ith()
	999	call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
[135]	1000	! Check:
	1001	! do i=1,iip1
	1002	! do j=1,jjp1
	1003	! do isoil=1,nsoilmx
	1004	! write(77,*) i,j,isoil," ",ith(i,j,isoil)
	1005	! enddo
	1006	! enddo
[1543]	1007	! enddo
[135]	1008
	1009
	1010
[1297]	1011	c slab_ocean_initialisation
	1012	c ------------------------------------------------
	1013	else if (trim(modif) .eq. 'slab_ocean_0') then
[1647]	1014
	1015	write(,) "No ocean yet on Titan ! Can't use this option"
	1016	stop
[135]	1017
[1297]	1018	else
	1019	write(,) ' Unknown (misspelled?) option!!!'
	1020	end if ! of if (trim(modif) .eq. '...') elseif ...
[135]	1021
	1022
	1023
	1024	enddo ! of do ! infinite loop on liste of changes
	1025
	1026	999 continue
	1027
[253]	1028
[1789]	1029	c================================================================
	1030	c Initialisation for methane surface tank (added by JVO 2017)
	1031	c================================================================
	1032	DO ig=1, ngridmx
	1033	tankCH4(ig)=2.0
	1034	ENDDO
	1035
[253]	1036	c=======================================================================
[135]	1037	c Correct pressure on the new grid (menu 0)
	1038	c=======================================================================
	1039
[588]	1040
[135]	1041	if ((choix_1.eq.0).and.(.not.flagps0)) then
	1042	r = 1000.*8.31/mugaz
	1043
	1044	do j=1,jjp1
	1045	do i=1,iip1
[588]	1046	ps(i,j) = ps(i,j) *
	1047	. exp((phisold_newgrid(i,j)-phis(i,j)) /
	1048	. (t(i,j,1) * r))
[253]	1049	end do
	1050	end do
	1051
[588]	1052	c periodicite de ps en longitude
[253]	1053	do j=1,jjp1
[588]	1054	ps(1,j) = ps(iip1,j)
[135]	1055	end do
	1056	end if
	1057
[588]	1058
[135]	1059	c=======================================================================
	1060	c=======================================================================
	1061
	1062	c=======================================================================
	1063	c Initialisation de la physique / ecriture de newstartfi :
	1064	c=======================================================================
	1065
[588]	1066
[135]	1067	CALL inifilr
	1068	CALL pression(ip1jmp1, ap, bp, ps, p3d)
[1889]	1069
[135]	1070
[1889]	1071	c=========================================================================
	1072	c Calcul de la dimension verticale pour la chimie - JVO 2017
	1073	c start_archive seulement, la grille verticale pouvant avoir ete modifiee
	1074	c==========================================================================
	1075
	1076	IF (choix_1.eq.0) THEN
	1077
	1078	! Calculate the # of upper chemistry layers with the "new" pressure grid
	1079	! For this we use Vervack profile for upper atmosphere with dz=10km
	1080
	1081	CALL gr_kim_vervack
	1082
	1083	WRITE(,)
	1084	WRITE(,) " With the compiled vertical grid we found :"
	1085	WRITE(,) " Number of upper chemistry layers =", nlaykim_up
	1086
	1087	! Regriding is then done, if needed
	1088
[1890]	1089	IF (callchim .and. nlaykimold.ne.nlaykim_up) THEN
[1889]	1090
	1091	WRITE(,) " Warning, nlaykimold=", nlaykimold
	1092	WRITE(,) ' which implies that a regriding on upper chemistry
[1890]	1093	& fields will be performed.'
[1889]	1094	WRITE(,)
	1095
[1890]	1096	! CALL regrid_kim(ngridmx)
[1889]	1097
	1098	ENDIF
	1099
	1100	endif ! of if (choix_1.eq.0)
	1101
	1102
[135]	1103	c-----------------------------------------------------------------------
	1104	c Initialisation pks:
	1105	c-----------------------------------------------------------------------
	1106
	1107	CALL exner_hyb(ip1jmp1, ps, p3d, beta, pks, pk, pkf)
	1108	! Calcul de la temperature potentielle teta
	1109
	1110	if (flagtset) then
	1111	DO l=1,llm
	1112	DO j=1,jjp1
	1113	DO i=1,iim
	1114	teta(i,j,l) = Tset(i,j,l) * cpp/pk(i,j,l)
	1115	ENDDO
	1116	teta (iip1,j,l)= teta (1,j,l)
	1117	ENDDO
	1118	ENDDO
	1119	else if (choix_1.eq.0) then
	1120	DO l=1,llm
	1121	DO j=1,jjp1
	1122	DO i=1,iim
	1123	teta(i,j,l) = t(i,j,l) * cpp/pk(i,j,l)
	1124	ENDDO
	1125	teta (iip1,j,l)= teta (1,j,l)
	1126	ENDDO
	1127	ENDDO
	1128	endif
	1129
	1130	C Calcul intermediaire
	1131	c
	1132	if (choix_1.eq.0) then
	1133	CALL massdair( p3d, masse )
	1134	c
	1135	print *,' ALPHAX ',alphax
	1136
	1137	DO l = 1, llm
	1138	DO i = 1, iim
	1139	xppn(i) = aire( i, 1 ) * masse( i , 1 , l )
	1140	xpps(i) = aire( i,jjp1 ) * masse( i , jjp1 , l )
	1141	ENDDO
	1142	xpn = SUM(xppn)/apoln
	1143	xps = SUM(xpps)/apols
	1144	DO i = 1, iip1
	1145	masse( i , 1 , l ) = xpn
	1146	masse( i , jjp1 , l ) = xps
	1147	ENDDO
	1148	ENDDO
	1149	endif
	1150	phis(iip1,:) = phis(1,:)
	1151
	1152	itau=0
	1153	if (choix_1.eq.0) then
	1154	day_ini=int(date)
	1155	endif
[1889]	1156
[135]	1157	c
	1158	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
	1159
	1160	CALL caldyn0( itau,ucov,vcov,teta,ps,masse,pk,phis ,
	1161	* phi,w, pbaru,pbarv,day_ini+time )
	1162
[588]	1163
[1415]	1164	CALL dynredem0("restart.nc",day_ini,phis)
	1165	CALL dynredem1("restart.nc",0.0,vcov,ucov,teta,q,masse,ps)
[135]	1166	C
	1167	C Ecriture etat initial physique
	1168	C
	1169
[1216]	1170	call physdem0("restartfi.nc",lonfi,latfi,nsoilmx,ngridmx,llm,
	1171	& nqtot,dtphys,real(day_ini),0.0,
	1172	& airefi,albfi,ithfi,zmea,zstd,zsig,zgam,zthe)
[1415]	1173	call physdem1("restartfi.nc",nsoilmx,ngridmx,llm,nqtot,
[1216]	1174	& dtphys,real(day_ini),
[1789]	1175	& tsurf,tsoil,emis,q2,qsurf,tankCH4)
[253]	1176
[135]	1177	c=======================================================================
[1543]	1178	c Formats
[135]	1179	c=======================================================================
	1180
	1181	1 FORMAT(//10x,'la valeur de im =',i4,2x,'lue sur le fichier de dema
	1182	*rrage est differente de la valeur parametree iim =',i4//)
	1183	2 FORMAT(//10x,'la valeur de jm =',i4,2x,'lue sur le fichier de dema
	1184	*rrage est differente de la valeur parametree jjm =',i4//)
	1185	3 FORMAT(//10x,'la valeur de lllm =',i4,2x,'lue sur le fichier demar
	1186	*rage est differente de la valeur parametree llm =',i4//)
	1187
[1216]	1188	write(,) "newstart: All is well that ends well."
	1189
[135]	1190	end
	1191

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

Download in other formats:

Original Format