Context Navigation

newstart.F @ 1647

Last change on this file since 1647 was 1647, checked in by jvatant, 8 years ago

+ Major clean of the new LMDZ.TITAN from too-generic options and routines (water, co2, ocean, surface type ...)
+ From this revision LMDZ.TITAN begins to be really separated from LMDZ.GENERIC
+ Partial desactivation of aerosols, only the dummy case is still enabled to keep the code running ( new aerosol routines to come in followings commits )

JVO

File size: 46.6 KB

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

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

Download in other formats:

Original Format